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Abstract 

We consider imbalanced Fermi gases with strong attractive interactions, for 
which Cooper-pair formation plays an important role. The two-component mix- 
tures consist either of identical fcrmionic atoms in two different hyperfine states, 
or of two different atomic species both occupying only a single hyperfine state. 
In both cases, the number of atoms for each component is allowed to be differ- 
ent, which leads to a spin imbalance, or spin polarization. Two different atomic 
species also lead to a mass imbalance. Imbalanced Fermi gases are relevant to 
condensed-matter physics, nuclear physics and astroparticle physics. They have 
been studied intensively in recent years, following their experimental realization 
in ultracold atomic Fermi gases. The experimental control in such a system 
allows for a systematic study of the equation of state and the phase diagram as 
a function of temperature, spin polarization and interaction strength. In this re- 
view, we discuss the progress in understanding strongly-interacting imbalanced 
Fermi gases, where a main goal is to describe the results of the highly controlled 
experiments. We start by discussing Feshbach resonances, after which we treat 
the imbalanced Fermi gas in mean-field theory to give an introduction to the 
relevant physics. We encounter several unusual supcrfluid phases, including 
phase-separation, gapless Sarma superfluidity, and supersolidity. To obtain a 
more quantitative description of the experiments, we review also more sophisti- 
cated techniques, such as diagrammatic methods and the renormalization-group 
theory. We end the review by discussing two theoretical approaches to treat the 
inhomogeneous imbalanced Fermi gas, namely the Landau-Ginzburg theory and 
the Bogoliubov-de Gennes approach. 
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1. Introduction 

Many fundamental discoveries have been made in the history of many-body 
quantum physics, which include new states of matter, unexpected phase tran- 
sitions, and novel macroscopic quantum effects. These discoveries have in com- 
mon that the rules of quantum mechanics are no longer restricted to the sub- 
atomic world, but also determine the collective behavior of systems that are 
large enough to be observed with the naked eye. The prime example that illus- 
trates such a discovery is the experimental observation of superconductivity by 
Kamerlingh Onnes [T]. In 1911, he found that the resistance of mercury sud- 
denly dropped to zero at a temperature of about 4 K. It was not immediately 
realized that a macroscopic quantum effect caused superconductivity, an impor- 
tant reason being that quantum mechanics was still at its infancy. Nearly half a 
century later superconductivity could be finally explained from first principles 
by Bardeen, Cooper and Schrieffer (BCS) [2J, who showed that in particular the 
pairing of electrons plays a crucial role. 

Pairing occurs in nature when two particles stick together due to an attrac- 
tive interaction. This attractive force may be strong, so that the pair is deeply 
bound and is hard to separate. The attraction may also be weak, so that the 
paired particles are much further apart, and a slight disruption can break the 
pair up. The group behavior of pairs can be quite different from their individ- 
ual behavior, which is particularly true for fermionic many-body systems. Here, 
pairing gives rise to a bosonic degree of freedom, whose statistics are freed from 
the Pauli principle that governs the behavior of the fermions. A macroscopic 
number of pairs can therefore occupy the same single-pair quantum state, lifting 
this state from the microscopic to the macroscopic world. The macroscopically 
occupied quantum state gives rise to the state of matter known as a Bose- 
Einstein condensate (BEC), which was already predicted by Einstein in 1925 
to occur in a noninteracting gas of bosons [3] . When Kapitsa, Allen and Mis- 
ener discovered flow without friction in the strongly interacting bosonic liquid 
4 He in 1938 |H[S], it took only a few months before a qualitative link between 
superfluidity and Bose-Einstein condensation was first made by London [5J. 

However, in fermionic systems it took some more time before the link be- 
tween superconductivity and condensation was noticed, since it was not until 
1956 that the microscopic BCS theory of superconductivity was finally formu- 
lated. This theory showed that a weak attraction between electrons, mediated 
by the subtle mechanism of lattice vibrations in metals (phonons) , caused a for- 
mation of loosely bound Cooper pairs, which could subsequently Bose-Einstein 
condense [2J. The BCS theory of pairing in Fermi mixtures revolutionized many- 
body quantum physics. It turned out to describe a very general mechanism that 
arises under a wide range of circumstances. An important reason for this is that 
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most matter in the world around us is made out of fermions, whether they are 
quarks, electrons, protons, neutrons, or more complicated composed particles, 
like fermionic atoms or molecules. Moreover, these fermions only need very 
little attraction to form a Cooper paired state. Indeed, Cooper showed that 
an infinitesimal amount of attraction is enough in the presence of a Fermi sea 
[7]. As a result, paired condensates have been observed or predicted to occur 
not only in a wide range of condensed-matter systems, but for example also in 
ultracold atomic gases El [10] or in the cores of neutron stars [UJ [12] . 

The BCS variational wavefunction itself also turned out to be much more 
general than initially expected. The many-body wavefunction was intended 
to describe the case of weak interactions between fermions that leads to the 
formation of loosely-bound pairs, whose size is much larger than the average 
interparticle distance. This regime is also called the BCS limit. As realized early 
on by Leggett [13], the BCS Ansatz also describes the case of strong attractions 
between the fermions leading to the formation of tightly bound molecules, whose 
size is much smaller than the interparticle distance. This regime is also called the 
BEC limit. Leggett predicted that the two extremes are continuously connected 
by a crossover, which was verified in recent years in a series of ground-breaking 
cold-atom experiments by the group of Jin using 40 K [9] and the groups of 
Ketterle, Thomas, Grimm, Salomon, and Hulet using 6 Li [HI \T5\ \T6 [ IT71 ITS] . 

The crossover experiments could be performed in the field of atomic quan- 
tum gases, because they have the unique feature that the effective interparticle 
interaction strength can be varied over an infinite range by means of a so-called 
Feshbach resonance (T9l [20] . In a Feshbach- resonant collision, two atoms col- 
lide and virtually form a long-lived molecule with a different spin configuration 
than the incoming atoms, after which the molecule ultimately decays into two 
atoms again. The two different spin configurations are also called the two differ- 
ent channels of the resonance. The scattering properties of the colliding atoms 
depend very sensitively on the energy difference of the molecular state with re- 
spect to the threshold of the two-atom continuum, which can be changed with 
an applied magnetic field. The theoretical prediction for the presence of tunable 
resonances in alkali atoms [21j and their actual realization in experiments [22] 
have enormously enhanced the power and scope of the experimental possibilities, 
especially in the case of ultracold Fermi mixtures. To illustrate this we notice 
that already shortly after Bose-Einstein condensation (BEC) was achieved in a 
dilute gas of bosonic alkali atoms by the group of Cornell and Wieman in 1995 
[2"3l , it was predicted that the superfluid regime could also be reached in a dilute 
mixture of fermionic atoms [8J. However, since the temperature for the BCS 
transition is exponentially suppressed for weak attractive interactions, the criti- 
cal temperatures initially turned out to be beyond the reach of experiments with 
ultracold Fermi mixtures. By using Feshbach resonances the interaction could 
be effectively enhanced which finally allowed for the observation of fermionic 
superfluidity and the BEC-BCS crossover. 

The crossover physics is usually studied in a two-component Fermi mixture 
with an equal number of atoms in each of the two different spin states. Then, 
a maximal number of Cooper pairs is created, which we can understand in 
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Figure 1: Illustration of the superfiuid-to-normal transition as a function of polarization 
P in a two-component Fermi mixture. In the balanced case, P = 0, attractive interactions 
between particles of different spin allow for the formation of a Cooper-pair condensate. Since 
the fully polarized case, P = 1, can be considered noninteracting, it leads to the normal state. 
Obtaining a detailed understanding of the system's evolution between these two limits is an 
important goal of this review. 



the following way. At low temperatures two particles predominantly interact 
through low-energy collisions for which the angular momentum is zero, also 
called s-wave collisions. Since for identical fermions the wavefunction should be 
anti-symmetric upon particle exchange, and since the two-particle wavefunction 
for s-wave scattering is symmetric in coordinate space, we conclude that s- 
wave interactions can only take place if the wavefunction is anti-symmetric in 
spin space. This is impossible to achieve if the particles only have access to 
one spin state. As a result, the fully polarized gas can typically be considered 
as noninteracting at very low temperatures, so that pairs between particles of 
the same spin are not formed. However, in the case of two spin states an anti- 
symmetric combination of the two spin states can be made and s-wave attractive 
interactions are possible. In the case of an equal amount of particles in both spin 
states, each particle can find another particle of a different spin to pair with, 
which leads to the formation of a supcrfluid after condensation of the pairs. 

From the above discussion it becomes immediately clear that something 
must happen as a function of spin imbalance or polarization, defined as P — 
(N + — iV_ ) / (N + + A_ ) with N a the number of particles in spin state a. Indeed, 
since at ultralow temperatures the gas is in a superfluid state for a spin-balanced 
mixture with attractive interactions, while it is noninteracting and thus in the 
normal state for a fully polarized mixture, it must undergo a phase transition as 
a function of polarization. This transition is illustrated in Fig. [I] The ultracold 
atomic Fermi gas was studied experimentally as a function of spin polarization 
for the first time by Zwierlein et al. and Partridge et al. [25], after which 
a large amount of experimental and theoretical activity followed. A main goal 
of this review is to achieve a detailed understanding of the phase diagram for 
the imbalanced Fermi mixture in the strongly interacting regime. Due to the 
generality of the pairing mechanism in fcrmionic many-particle systems, this 
topic is also of direct interest to condensed-matter, nuclear and astroparticle 
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physics [251 H3]. 

Although in the past century an impressive amount of knowledge about 
pairing in Fermi mixtures has been achieved, many questions still remain. An 
important example is the phenomenon of high-temperature superconductivity, 
which was discovered in a certain class of ceramic materials in 1986 [27|. Here, 
the microscopic nature of the pairing has turned out to be one of nature's best 
kept secrets. In order to address the many remaining fundamental open ques- 
tions that involve pairing in Fermi mixtures, it is extremely helpful to have a 
physical system that is experimentally very clean, highly tunable and easy to 
probe. Such a system then gives rise to a well-defined microscopic hamiltonian 
where the interactions between the particles are accurately known. It is thus 
an ideal starting point for calculations of thermodynamic functions that govern 
the macroscopic properties of the system. These calculations go hand in hand 
with detailed experimental studies, where either experiments set benchmarks 
for sophisticated theoretical studies, or where theoretical predictions give rise 
to landmark experimental discoveries. In this way, a much more detailed under- 
standing of fermionic superfluidity is achieved. Such a system is indeed nowa- 
days available, namely in the field of ultracold atomic quantum gases. Whether 
the mystery of high-temperature superconductivity will ultimately be unrav- 
elled in the world of cold atoms still remains to be seen, but many fundamental 
discoveries have already followed each other up at high speed, while many more 
arc likely to follow soon [H? l 125 1 121? 1 150] . 

1.1. Ultracold atomic quantum gases 

When Bosc-Einstein condensation (BEC) was ultimately realized in trapped 
gases of bosonic atoms [23J [3lJ [32] , a completely new category of systems for 
studying macroscopic quantum effects became available. A crucial ingredient 
for reaching BEC was the development of laser cooling, with which atoms are 
made to absorb and emit photons such that a momentum transfer takes place 
[551 1531 [551 [55] , This was used to create a friction force for atoms, slowing 
them down to a near standstill and thereby making the gas cloud enter the 
microKelvin regime. Equally important was the trapping of the atoms in a 
magnetic trap [57]. As a result, the gas could be held together without making 
contact to material walls that cannot be cooled to such ultralow temperatures. 
Moreover, the depth of the magnetic trap was easily lowered which allows the 
most energetic atoms to escape. If the remaining atoms undergo elastic colli- 
sions, then they re-thermalize at a lower temperature, allowing for an evapo- 
rative cooling of the trapped quantum gas [38 [ 36] . This mechanism is similar 
to the way in which a cup of coffee gets cold during lunch time. Evaporative 
cooling turned out to be extremely efficient, allowing experimentalists to reach 
a temperature of 170 nK at a particle density of 2.5 x 10 12 cm' 3 , where BEC 
was finally observed [23]. Note that this means that an ultracold atomic quan- 
tum gas is typically about ten million times thinner than air. Moreover, it also 
means that a condensed gas of alkali atoms is strictly speaking in a metastable 
state, because at these ultralow temperatures the stable state is actually a solid. 
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However, the solid formation takes so long on the time scales of the experiment, 
that we can consider the dilute gas to be in thermodynamic equilibrium. 

Although the discovery of the Bose-Einstein condensed state of matter was 
extremely exciting in itself, it turned out to be only the beginning for the explo- 
ration of new macroscopic quantum effects in ultracold gases. A reason for the 
fruitfulness of atomic many-particle systems is their extreme cleanness, mean- 
ing that there are essentially no impurities in the experiments. This is in a 
sharp contrast to solid-state systems, where the understanding of many-body 
experiments can be extremely severed due to an uncontrolled amount of imper- 
fections in the materials. Disorder in atomic gases is typically negligibly small, 
unless it is deliberately added [3JE HOI SI] , which led to a very controlled study 
of Anderson localization. 

This brings us automatically to a second reason for the fruitfulness of degen- 
erate quantum gases, namely the amazing amount of experimental control that 
is currently achieved [TUl HSj ■ We have already mentioned the two-channel Fes- 
hbach resonance with which the interatomic interaction strength can be tuned. 
Another important example of control is the external trapping potential of the 
gas cloud, which can be precisely tailored experimentally. A particularly conve- 
nient setup is achieved when the atoms are optically trapped with the use of the 
strong electric fields in laser beams. The laser beams can then be made counter- 
propagating which leads to an intense standing wave of light. This creates a 
periodic potential for the atoms due to the Stark effect, giving rise to a so-called 
optical lattice [32 S3] - Optical lattices can be used to simulate ionic lattices, 
which offers the opportunity to explore various aspects of condensed-matter 
physics in the very clean environment of ultracold atoms |28j . The depth of the 
periodic potential is now easily tunable by varying the laser intensity, while the 
period of the lattice is directly related to the frequency of the laser light. An 
even more flexible possibility to create arbitrary potential landscapes for the 
atoms is achieved by shining a laser onto a holographic mask [44 . . 




In the so-called Hubbard model [45] , a paradigm in condensed- matter physics, 
particles are allowed to tunnel between adjacent lattice sites and they have an 
on-site interaction. When an interacting atomic gas at ultralow temperatures 
is loaded into an optical lattice, then an essentially perfect manifestation of the 
Hubbard model is realized. This was first realized by Jaksch et al. [15], who 
treated repulsively interacting bosons. Considering the case of an equal number 
of atoms and lattice sites, we have that the bosons are expected to form a de- 
localized superfluid state when the tunneling strength t is dominant. However, 
when the on-site interaction strength U is dominant, the ground state is given by 
the so-called Mott-insulator state, which has precisely one localized atom at each 
lattice site. At zero temperature, a transition occurs between these two phases 
as a function of t/U , which is solely driven by quantum fluctuations and not 
by thermal fluctuations [47] ■ This quantum phase transition was observed in a 
landmark atomic physics experiment by Greiner et al. |43| . Here, the tunneling 
strength was varied by simply changing the laser intensity. Note that an analo- 
gous control over the tunneling would be hard to achieve in a condensed-matter 
system. An ultracold Fermi mixture can also be loaded into an optical lattice, 
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Figure 2: Impression of several configurations that are realizable with optical lattices, a) If 
an intense standing wave of laser light is applied in one direction, then the potential landscape 
for the atoms becomes a stack of quasi-two-dimensional 'pancakes', b) If a second standing 
wave is added in a perpendicular direction, then an array of quasi-one-dimensional 'cigars' 
arises, c) If counter-propagating laser beams are added in the third direction, then a lattice 
of point-like sites arises. 



which leads to a realization of the fermionic Hubbard model. This has led to 
the observation of the Mott-insulator phase with ultracold fermions [48,09]. In 
the absence of doping the ground state of the Hubbard model is predicted to 
be the anti-ferromagnetic Neel state, while the doped fermionic Hubbard model 
has been conjectured to encompass high-temperature d-wave superconductivity 
|50j . It would be exciting to study these phenomena in the controlled environ- 
ment of ultracold atomic gases. However, to this end a further decrease in the 
temperature must be achieved experimentally, which is a very challenging task. 

There are many more interesting applications of optical lattices. For ex- 
ample, they can be used to create low-dimensional quantum gases. Consider 
a lattice which is very steep in, let's say, the x direction, such that the ultra- 
cold particles can only occupy the lowest-lying quantum state with an energy 
Huj x /2 in this direction. Then, the kinetic degrees of freedom are frozen out 
in the x direction, which effectively lowers the dimension of the system. In 
particular, with a very steep optical lattice in two directions it is possible to 
create a two-dimensional array of effectively one-dimensional tubes, while a 
very steep optical lattice in one direction leads to a one-dimensional stack of 
two-dimensional pancakes, as shown in Fig. [2] Low-dimensional quantum gases 
give rise to rise to intriguing strongly-correlated behavior that is very different 
from the three-dimensional case |51j . Experimentally observed examples include 
the exotic Bcrczinskii-Kosterlitz-Thouless phase transition in two dimensional 
gases [52] [53] [54] , the coherence dynamics of one-dimensional Bose gases [55] 
and exotic pairing of spin-imbalanced Fermi gases in one dimension [56 . 

But the control in ultracold atomic gases goes further. Feshbach resonances 
are for example not restricted to s-wave scattering, so that also resonances of 
higher angular momentum have been observed. This has led to the formation 
of p-wave Feshbach molecules [571 EH EH] , which upon condensation would give 
rise to superfluidity in a nonzero angular momentum state. Analogous to liq- 
uid 3 He, this supcrfluid has a more complex order-parameter structure |60l 161] , 
which may give rise to exotic features such as the presence of Majorana fermions 
in vortex excitations, topological phase transitions and topologically protected 
quantum computing [52J [5U [S3] . A different research direction is the creation 
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of condensates, where the particles interact through long-ranged anisotropic 
dipole-dipole interactions. Dipolar effects have already been observed in con- 
densates of 52 Cr atoms [64], whose magnetic moment leads to weak dipole in- 
teractions. The interaction is much stronger for heteronuclear molecules that 
have a permanent electric dipole moment. An ultracold gas of ground-state 
40 K 87 Rb molecules has been achieved by associating 40 K and 87 Rb atoms us- 
ing a Feshbach resonance and then optically pumping the 40 K 87 Rb dimers to 
their rovibrational ground state [65]. The interesting possibilities with dipolar 
molecules vary all the way from the stabilization of supersolid phases 66 to the 
possibility of quantum computing |67j . Another promising line of research is 
to create effective artificial vector potentials for cold atoms using Raman tran- 
sitions [55], which have among others resulted in artificial magnetic fields [55] 
and spin-orbit coupling ;70 for ultracold neutral atoms. As a result, the rich 
physics of gauge theories for charged particles can be studied in the versatile 
atomic physics environment. A final example of control that is of particular 
relevance to this review is the controlled preparation of an atomic gas in a se- 
lective set of internal quantum states. In a two-component 6 Li mixture it is 
possible to precisely control the atom number in each of lithium's two lowest 
hyperfine states by inducing nuclear spin flips [211 [25]. As a result, the fun- 
damental phase diagram of the two-component Fermi mixture can be studied 
as a function of temperature, polarization and interaction strength, of which 
the three-dimensional case is a central topic of this review. Moreover, if two 
different atomic species are used, then a fourth axis enters the phase diagram, 
namely that of a mass imbalance. 

To make a long story short, the possibilities seem endless and the field of ul- 
tracold atoms is able to address fundamental questions about many-body quan- 
tum physics in great detail. Therefore, atomic quantum gases are sometimes 
also referred to as ideal quantum simulators. They allow for systematic stud- 
ies of an enormous variety of hamiltonians, ranging from weakly interacting to 
strongly interacting, from one dimensional to three dimensional, from disordered 
to clean, from homogeneous to periodic, where the microscopic parameters are 
always precisely known and widely tunable. 

1.2. This review 

In this review, we consider three-dimensional two-component Fermi gases 
with attractive interactions so that pairing is expected to occur. However, 
the equal-density Cooper pairing mechanism is frustrated when an increasing 
spin imbalance is introduced to the system. We are interested in the way the 
system responds to this kind of frustration, where multiple scenarios might 
be realized, as we discuss in the next section. In recent years, many important 
experiments have been performed on the imbalanced Fermi gas, after this system 
was first realized experimentally in the field of ultracold quantum gases at MIT 
by Zwierlein et al. [H] and at Rice University by Partridge et al. 25 . The 
initial studies of these groups focussed on experimentally determining the phase 
diagram of the strongly interacting Fermi g clS £IS cl function of spin polarization 
and temperature [71] [72]. Later, the imbalanced atomic Fermi gas was also 
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experimentally studied at ENS, where Nascimbene et al. accurately measured 
the collective modes 73 and the equation of state of the imbalanced Fermi gas 
[111 [75]. In most recent experimental studies, also nonequilibrium aspects of 
the imbalanced Fermi gas are considered, such as spin transport [751 IZZ] and 
metastable states [75] . 

An early excellent review that focusses on the properties of imbalanced Fermi 
gases in mean-field theory at zero temperature along the BEC-BCS crossover 
can be found in Ref. |79j . Other shorter, but very insightful reviews on the 
imbalanced Fermi gas are Refs. 80, 81] . In the present long review, we mainly 
focus on the equilibrium properties of imbalanced Fermi gases in the unitarity 
limit. A central difference with Ref. [79] is that we spend much attention to 
the effects of interactions beyond mean-field theory by discussing diagrammatic 
methods and renormalization-group theory. Moreover, we consider the effects 
of the trap beyond the local-density approximation by discussing the Landau- 
Ginzburg approach and the Bogoliubov-de Gennes equations. Finally, we also 
treat nonzero temperatures and the mass-imbalanced case. The phase diagram 
of the imbalanced Fermi gas is a topic of fundamental interest to various areas 
of physics |26L I12j . where we can think of neutron-proton pair condensation 
in nuclear physics [82] . or color superconductivity of strongly degenerate quark 
matter in the core of neutron stars |11] . We start out by qualitatively discussing 
which kind of paired phases might occur. 

1.2.1. Unusual superfluid phases 

The 'standard' way in which fermions form a superfluid is when identical 
particles (such as electrons) of different spin undergo s-wave attractive inter- 
actions and form pairs [2]. These Cooper pairs constitute a bosonic degree of 
freedom, and at low enough temperatures they condense forming a superfluid. 
We note here that although in three dimensions the formation of a Bose-Einstein 
condensate typically goes hand in hand with the onset of superfluidity, they are 
not identical phenomena. For example, the enhancement of fluctuation effects 
in lower dimensions can destroy a condensate while keeping the property of 
superfluidity intact. Moreover, in three dimensions the ideal Bose gas is fully 
condensed at zero temperature, but it is usually not considered to be a superfluid 
[S3]. The latter is because its critical velocity, which is the highest velocity for 
the gas to flow without friction, is then equal to zero. In this review we consider 
the interacting three-dimensional case, so that we use the terms Cooper-pair 
condensate and superfluid state interchangeably. 

When there is a population imbalance between the two spin states of the 
Fermi mixture, not all particles can find a partner to pair up with. In this re- 
view we will study the phase diagram of the spin-imbalanced mixture for strong 
attractive interactions, both in the presence and the absence of a trapping po- 
tential. Due to the enhanced richness of this system, we should consider at least 
two additional superfluid phases that are not present in the spin-balanced case. 
First of all, we find the occurrence of a phase-separated phase at temperatures 
close to zero. This phase can arise because there is a discontinuous or first-order 
transition in the system between a superfluid state that has a small polarization 
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and a normal state that has a high polarization [75] . The phase-separated 
state then occurs when the total polarization of the system is in between these 
two extremes. As a result, the system spatially separates into a supcrfluid do- 
main and a normal domain, between which there is a first-order interface, i.e. a 
domain wall. In the presence of a trap, the superfluid domain is then found in 
the center of the trap, while the normal domain surrounds it. 

The second less-standard superfluid phase that we consider is the so-called 
Sarma phase (SU [85] , which has a gapless excitation spectrum. As a result, 
this superfluid phase is polarized even at zero temperature. It is also called the 
interior-gap phase, or the breached-pair phase [BE1 IS3 [55] • The Sarma phase is a 
supcrfluid with one or more Fermi surfaces, so that it is different from the gapped 
BCS phase which has no Fermi surface, and the two phases are separated by 
a quantum phase transition at zero temperature. For the mass-balanced case, 
we see later on that mean-field theory predicts the Sarma phase to become 
unstable at very low temperatures, preventing the observation of the quantum 
phase transition. At nonzero temperatures, the distinction between the Sarma 
phase and the BCS phase is not sharp because temperature causes both phases 
to be polarized and to not have a sharp Fermi surface. As a result, the gapless 
regime and the gapped regime are then connected by a crossover |89] . In the case 
of a large mass-imbalance, the Sarma phase becomes stable at zero temperature 
according to mean- field theory [55] • For completeness, we mention also other 
exotic paired phases that have been predicted to enter the phase diagram of the 
spin imbalanced Fermi gas, namely states with a deformed Fermi surface [SO] , 
p-wave superfluidity [91] [92] , and supersolid phases [93] [94] [80] . 

Although in this review we primarily study the Fermi gas in the unitarity 
limit as a function of temperature and spin polarization, there are more param- 
eters to be varied. The phase diagram can for example also be explored as a 
function of interaction strength, and as a function of the mass imbalance between 
the two different components of the mixture. The situation of a mass imbalance 
is currently being studied by many experimental groups, where in particular 
the 6 Li- 40 K mixture seems promising |95l 1961 1971 198| . However, pair condensa- 
tion with different fermion masses has not been observed yet. The mean- field 
phase diagram of the 6 Li- K mixture not only encompasses the rich physics 
from the mass-balanced case, but it has even more structure [§S] • Namely, next 
to the presence of Sarma physics and phase separation, mean-field theory also 
predicts a Lifshitz point in the phase diagram for the unitarity limit, as we see 
explicitly in Section |3.7| A Lifshitz point also enters the phase diagram in the 
mass-balanced case, but then at weak interaction strengths, making the corre- 
sponding Lifshitz temperature very low |100l 110 Jl I102j . At the Lifshitz point, 
the effective mass of the Cooper pairs turns negative, which signals a transition 
to an inhomogeneous superfluid. This exotic possibility was first investigated 
by Larkin and Ovchinnikov (LO) [103] . who considered a superfluid with a 
standing-wave order parameter, namely A(R) = AqCos(K • R). The order pa- 
rameter spontaneously breaks translational symmetry, and because the atom 
densities depend on the absolute value of the superfluid order parameter, the 
densities also oscillate in space. Moreover, the phase is typically energetically 
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more favorable than a plane-wave order parameter A(R) = Aoe iK R , which was 
considered independently by Fulde and Ferrell (FF) (104] . The FF phase leads 
to a propagating superfluid that spontaneously beaks time-reversal symmetry, 
while it does not give rise to oscillating particle densities. 

The combination of superfluidity with crystalline order in the densities makes 
the LO phase a particular supersolid |105llM] . To be more precise, the LO phase 
can also be called a density wave of a Cooper-pair condensate [1061 11071 11081 
109 , since the order parameter is a product of a pairing amplitude Ao and a 
standing wave cos(K ■ R). This is in contrast to the supersolid phases that are 
usually considered, for which the superfluid order and the crystalline order are 
completely independent of each other. A special property of the pair-density 
wave is that the order parameter structure not only allows for the occurrence of 
integer- valued vortices (as for ordinary superfluids) and the occurrence of dislo- 
cations (as for ordinary density waves), but also for composite objects that are 
half vortex and half dislocation [10711109] . A careful study of the low-energy fluc- 
tuations of the LO phase, which breaks continuous rotational and translational 
symmetry, reveals that it is actually a phase with only algebraic long-range or- 
der, also known as a superfluid smectic liquid crystal [109] . Moreover, upon 
melting of the crystal other exotic phases can arise, such as fractionalized Fermi 
liquids and superfluids of pairs of Cooper pairs |109] . The physics of the FF and 
LO phases have intrigued the condensed-matter community for many decades, 
but it has been experimentally very challenging to prove the existence of these 
phases unambiguously [13] . The same holds for supersolidity in general, where a 
famous example is the heavily debated experiment with 4 He by Kim and Chan 
[110] , A promising way to observe FFLO phases directly is to study the one- 
dimensional spin-imbalanced Fermi gas [56] . for which the FFLO phases have 
been predicted to occupy large parts of the phase diagram [1111 11121 1113] . 

1.2.2. Outline 

Having mentioned several superfluid phases of interest that we encouter in 
this review, we now give a more precise outline of the topics that are to be 
treated. We begin in Section [2] with a comprehensive discussion of two-body 
scattering in ultracold atomic gases. This is necessary for understanding in more 
detail the interaction mechanisms on the two-body level before trying to tackle 
the many-body properties. In particular, we focus in this section on Feshbach 
resonances, which are so crucial for ultracold atom experiments these days. 
We end the section with a qualitative discussion of an important experimental 
example that reveals the power of using Feshbach resonances, namely the study 
of the BEC-BCS crossover. Having discussed the two-body physics, we turn 
to the central topic of this review in Section [3j namely the imbalanced Fermi 
gas. We start out by giving a historical overview of the early atomic-physics 
experiments that have been performed with this system up to the present status 
of the field. This overview will give us a taste of the physical phenomena that we 
wish to explain. Next, we discuss extensively the simplest theoretical approach 
to describe pairing in Fermi systems, namely the Bardeen-Cooper-Schrieffer 
(BCS) mean-field theory. We find that this simple approach can already account 
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qualitatively for many features observed in the experiments done in the strongly- 
interacting regime. We calculate the mean-field phase diagram for the solely 
spin-imbalanced case in the unitarity limit both for the homogeneous case and 
for the trapped case. For the latter, the local-density approximation is applied. 
We end this section with the mean-field phase diagram of the mass-imbalanced 
6 Li- 40 K mixture. 

Although mean-field theory is useful for obtaining a qualitative understand- 
ing of the relevant physics in imbalanced Fermi gases, it is not very accurate 
quantitatively in the unitarity regime. In Section [4] we discuss diagrammatic 
approaches that go beyond mean-field theory and that considerably improve 
agreement with experiments and Monte-Carlo calculations. First of all, we 
treat the Gaussian fluctuations of the BCS order parameter, also known as the 
Nozieres-Schmitt-Rink (NSR) approximation, which can be evaluated exactly 
and gives a significant improvement in the description of the BEC-BCS crossover 
for balanced Fermi gases compared to mean-field theory. For example, it leads 
to a rather accurate description of the critical temperature curve as a function 
of interaction strength, whereas mean-field theory fails to describe the critical 
temperature on the BEC side. Despite this remarkable success, the NSR method 
unfortunately fails upon application to the imbalanced Fermi gas. Therefore, 
we also discuss other techniques. First of all, we calculate the self-energy of 
fermions using a many-body transition-matrix approach. This self-energy can 
be consequently incorporated into an equation of state that is in very good agree- 
ment with experiments and Monte-Carlo simulations. Moreover, in Section [5] 
we apply the wilsonian renormalization scheme to the imbalanced Fermi gas in 
two different ways. Namely, first we apply it directly to the fermionic action, 
in order to not only include the effect of fermionic selfenergies, but also of the 
screening of the interaction by particle-hole excitations. We use this method to 
calculate more accurately the tricritical point of the phase diagram in the uni- 
tarity limit. We also apply the wilsonian renormalization scheme to the action 
for the Cooper pairs, which can be exactly obtained from the fermionic action 
by a so-called Hubbard-Stratonovich transformation. We find that by including 
the interaction between the Cooper pairs, the problems encountered with the 
Nozieres-Schmitt-Rink approximation are solved. 

Having obtained a rather detailed understanding of the physics beyond 
mean-field theory for the homogeneous imbalanced Fermi gas, we turn to the in- 
homogeneous case. In Section [6] we describe the trapped Fermi gas beyond the 
local-density approximation (LDA), because this approximation breaks down 
near the superfluid-normal interface that is observed for low-temperatures in 
experiments. We describe two approaches beyond LDA, namely the Landau- 
Ginzburg theory and the Bogoliubov-de Gennes equations. We use the Landau- 
Ginzburg approach to compare with experiments and to calculate the resulting 
surface tension. When using the Bogoliubov-de Gennes approach a proximity 
effect is observed, meaning that there are oscillations in the order parameter 
near the superfluid-normal interface, which are absent in the Landau-Ginzurg 
approach. Finally, we end Section[7]to summarize our conclusions and to discuss 
some issues that are still open and that require further research. 
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This review is meant to be introductory and pedagogical in style, meaning 
that it is aimed to be comprehensible by graduate students (both theoretical and 
experimental) entering this particular field of research. Only Sections [4] and [5] 
request previous knowledge of quantum field theory in the functional formalism 
and are therefore mainly focussed on the theoretically interested reader. The 
other parts can easily be followed without reading these two sections in great 
detail. Whenever steps in derivations are skipped we have tried to include a 
reference where the corresponding steps are explicitly taken. We hope that this 
approach has resulted in a useful review for the reader. 

2. Scattering and Feshbach resonances 

We start our discussion of interacting atomic Fermi gases by describing in 
some detail the relevant two-body physics of these systems, which is necessary 
as an input for the many-body theory. To this end we introduce several useful 
concepts from general two-body scattering theory. The machinery to deal with 
scattering problems is well established and can be found in most textbooks on 
quantum mechanics, see e.g. Refs. |114l 11151 1116j . Here, we only briefly state 
the results that are of most relevance to us. The main goal of this section is 
to understand the physical properties of Feshbach resonances, which have been 
crucial to all experiments with imbalanced Fermi gases. More extensive reviews 
on Feshbach resonances can be found in Refs. |117|I118] , We end the section by 
briefly describing qualitatively an important application of Feshbach resonances, 
namely the study of the BEC-BCS crossover. 

2.1. Single- channel scattering 

Consider two particles of mass m that elastically scatter from each other 
under the influence of an interaction potential V(r) that depends on the relative 
coordinate r of the particles. If the Schrodinger equation is separated into a part 
describing the center-of-mass motion and a part describing the relative motion 
of the two particles, then the center-of-mass part behaves like a single free 
particle of mass 2m, while the relative part behaves like a single particle with 
reduced mass m/2 moving in the potential V(r). Starting point is the relative 
Schrodinger equation 



where Hq = — /i 2 V 2 /m for the relative kinetic-energy operator. We start with 
looking at solutions where the particles enter the scattering region in a plane 
wave state |k) with energy E = H 2 k 2 /m = 2ek, which is conserved in the elastic 
collision. The state |k) also solves the Schrodinger equation in the absence of 
the interaction potential. Since Eq. ([lj is the time-independent Schrodinger 
equation, it may be interpreted as describing a steady-state solution for a con- 
tinuous stream of particles in state |k) scattering from a potential. It turns 
out to also describe the scattering of wavepackets, provided that the typical 
wavelength of the packets is much larger than the range of the interaction [115] . 




(1) 
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This condition is typically well satisfied for dilute ultracold atomic gases, which 
are characterized by large de Broglie wavelengths and by interactions that are 
very short ranged compared to the average interparticle distance. 

Eq. (lj) can be formally solved by introducing scattering states |V>k + ^) that 
satisfy the following recursion relation, also known as the Lippmann-Schwinger 
equation, 

IVk / — l"7 "T" ^o\**k.j y ivk 

where Go(2ek) = (2ek — Ho + «0) _1 with the infinitesimal positive imaginary 
part being the standard way to treat the singular nature of this operator. In 
a time-dependent formulation, this latter procedure follows immediately from 
demanding that the particles were free in the remote past |115j . In scatter- 
ing theory a central role is played by the two-body transition matrix, defined 



|^») = |k) + G„(2 £k )^[ +) ), (2) 



through F|V k +) ) = T 2b (2e k )|k). As a result, Eq. ([2) leads to 

f 2b (2e k )|k) = {V + yGo(2 ek )T 2b (2e k )} |k). (3) 

This form of the Lippmann-Schwinger equation can also be generalized to an 
operator equation, whose recursive solution gives rise to the so-called Born series 

f 2h {E) = V + VG Q (E)V + VG (E)VG (E)V +..., (4) 

for a certain energy E. More shortly, 

f 2b {E) = V + VG(E)V, (5) 

where G(E) = (E — H — V)^ 1 . The last equation indicates that singularities 
in the transition matrix are found at the exact eigenenergies of the two-body 
problem. 

More insight is gained by realizing that away from the scattering region 
the influence of the scattering potential becomes negligible, which leads to a 
noninteracting radial Schrodinger equation. Here, we demand that the solution 
describes an incoming plane wave with wavevector k and a freely propagating 
outward radial flow of probability flux. It is possible to show that the flux- 
normalized solution at large r is given by, see e.g. |115j . 

^ k +) (r)^ k ' r + / 2b (k',k)^, (6) 
where k' = kr . The two-body scattering amplitudes then satisfy 

/ 2b (k',k) = - 4 ^(k^k' k +) > = - 4 ^(k'|T 2b (2 ek )|k) 



OC 

e=o 



T 2b (fc)iMcostf). (7) 
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In the second line the partial-wave method was employed to expand the scat- 
tering amplitude and the transition matrix, where Pi{x) are the Legendre poly- 
nomials and ■& is the angle between k and k'. We note that in obtaining the 
expression of Eq. ^ the normalization (r|k) = e lkr was used for the plane 
waves. This reveals that the ket |k) itself is not dimensionless, which is a con- 
sequence of using the continuum limit for the states in k space. Moreover, the 
expansion in Eq. ^ is useful for a spherically symmetric interaction, for which 
the scattering problem is symmetric around the axis of incidence. From the 
conservation of angular momentum and probability flux, it can be shown that 

f " h{k) = kcot[5 e (k)]~tk' (8) 

where Si(k) is the phase shift of the £-th partial wave due to the elastic scatter- 
ing. 

At low energies, the phase shift is governed by the expansion [114] 

k 2e+1 cot[5 e {k)] = - (-) + 0(k 2 ), (9) 

so that for small wavevectors the partial scattering amplitudes are given by 
ff°{k) ~ — a 2i+1 k 2e with at the scattering parameter of the ^-th partial wave. 
For £ = 0, ciq is called the s-wave scattering length, which from now on we sim- 
ply denote as a, since we only consider scattering at zero angular momentum. 
This also avoids confusion with the same symbol often used for the Bohr radius. 
From the behavior of the partial scattering amplitudes, we thus see that at low 
momenta the s-wave scattering scattering with zero angular momentum (£ = 0) 
is dominant. The reason why collisions with higher angular momentum are sup- 
pressed is also readily understood from a more physical point of view. Namely, 
for £ 7^ the relative radial Schrbdingcr equation contains a repulsive centrifugal 
term h £(£ + l)/mr 2 , which acts as an energy barrier. At low temperatures the 
particles do not have enough kinetic energy to overcome this barrier and thus 
the higher-order partial waves do not feel the short-ranged interaction. These 
waves behave as if the interaction potential were not there, giving rise to quickly 
vanishing phase shifts. For s-wave scattering, we have that the scattering length 
is given by 

a= _ lim w (10) 

fc->o k 

Moreover, we find from Eq. {8} that / (0) = -a and from Eq. Q that T 2b (0) = 
ATrh 2 a/m. In the specific case of a hard-core potential, the scattering problem 
can easily be solved analytically, which leads to the s-wave scattering length 



being equal to the radius of the impenetrable core |117j . Therefore, a positive 
s-wave scattering length can be interpreted as the the effective hard-core radius 
of the interaction. 

2.2. Two-channel scattering 

We thus found that at low energies the two-body transition matrix is dom- 
inated by the s-wave scattering length, which is interpretable as giving the 
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effective interaction strength. Therefore, the scattering length is an important 
quantity to know in the description of ultracold atomic gases. An accurate 
calculation of the scattering length from first principles can be a challenging 
task, because typically the atomic interaction potentials are not known pre- 
cisely enough. Fortunately, the scattering length can be directly measured by 
experiments. This was for example done by Inouye et al. [55], who showed 
experimentally that the effective interaction strength can be tuned over a large 
range by varying the applied magnetic field. The mechanism that allows for 
this extremely useful control is nowadays called a Feshbach resonance [19] . and 
recent extensive reviews about this topic can be found in Rcfs. [1171 1118] , In 
the following discussion we mainly follow the lines of Ref . |117j . 

In an atomic Feshbach-resonant collision, the two incoming atoms virtually 
form a molecule with a different spin configuration, after which this molecule 
decays into two atoms again. The scattering properties depend very sensitively 
on the energy difference between the molecular state and the threshold of the 
two-atom continuum. This energy difference can be controlled with an applied 
magnetic field, because the difference in spin states between the incoming atoms 
and the molecule gives rise to a different Zeeman shift. To further study the 
physics, we use a two-channel model as illustrated in Fig.[3j The atoms approach 
each other in the atomic channel, while there is a bound state close to the 
threshold of the atomic continuum in the molecular channel. We are interested 
in the regime where the collisional energy is lower than the Zeeman splitting. 
Due to energy conservation the molecular channel is then inaccessible for the 
atoms at large separations from each other after the collision. Therefore, the 
molecular channel is also called the closed channel, while the atomic channel 
is called the open channel. We consider only the bound state in the molecular 
channel that is closest to the atomic continuum, since this state is dominant in 
affecting the scattering process. 

In the absence of a coupling between the channels, the molecule is described 
by the wavefunction \ip m ), which we call the bare molecular state, and has an 
energy 5^, which we call the bare detuning. The bare molecular state and the 
bare detuning are only a solution to the uncoupled Schrodinger equation in the 
closed channel. This is to be contrasted with the so-called dressed molecular 
state, which is the bound-state solution when the two channels are coupled. In 
the open channel the two atoms interact through a short-ranged potential Vbg, 
which depends on their relative coordinate r and which is called the background 
interaction. The atomic and the molecular channel are coupled by V am due to the 
possibility of electronic spin flips caused by hyperfine interactions. To find the 
energy and the wavefunction of the dressed molecule, the two-body problem is 
separated into a center-of-mass part and a relative part. The center-of-mass part 
gives rise to a trivially solvable free-body problem, while the relative Schrodinger 
equation becomes 

Ho + Vbg V am 

Vam 5 h 

where \^p a ) denotes the wavefunction in the atomic channel, IV'm) is the wave- 



Z m \4> m ) 
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y/1 - Z m |V>a) 
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Figure 3: Illustration of the two-channel structure of a Feshbach resonance. Shown are the 
atomic or open channel and the molecular or closed channel in the absence of a coupling 
between the two. The interaction potential in the open channel is called the background 
interaction Vb g (r), which depends on the interatomic separation r. The relevant feature of 
the closed channel is a bound state, also called the bare molecular state, that has a small 
energy difference 6^ with the atomic continuum. Due to a different spin configuration of the 
atoms in the two channels, there is a Zccman shift of A/iB with Afi the difference in magnetic 
moments and B the magnetic-field strength. 



function of the bare molecular state, Z m is the so-called wavefunction renor- 
malization factor that normalizes the wavefunction of the dressed molecule to 
unity, and Hq = —h 2 'V^/'m is again the relative kinetic energy operator. 

In this section we consider only s-wave scattering, which is dominant at 
low energies. As explained before, s-wave scattering of fermions only occurs 
between particles in two different internal states. In the experiments with 6 Li of 
interest to us, a two-component mixture is prepared in the lowest two hyperfinc 
states of the ground state. They are denoted by |1) and |2), which are shown 
in Fig. |4^b). Ultracold 6 Li is an experimentalist's favorite due to the easily 
accessible Feshbach resonance at a magnetic-field strength of Bq = 834 G, which 
is an extremely broad resonance. The resonance is shown in Fig. |4]Ja). Near 
834 G, the two lowest lying hyperfine states have the single valence electron of 
lithium anti-aligned with the applied magnetic field, i.e. m s = —1/2 [H31 1119] . 
These two states then differ in their projection of the nuclear spin. If two 
incoming atoms in state |1) and \2) collide, then the electronic spins of the 
atoms point in the same direction, so that the open channel gives rise to a spin 
triplet potential. During the collision an electronic spin can be flipped, which 
brings the two atoms in a closed spin singlet channel. The difference in magnetic 
moments between the two channels is A/j, — 2fi^ with /iB the Bohr magneton. 
As we increase the magnetic-field strength B from the resonance strength Bq, 
then the triplet potential is lowered in energy by an amount 8 = A/i(B — Bq), 
which shifts the atomic continuum down with respect to the bare molecular 
bound state as shown in Fig. [5] Equivalently, we can also say that the bare 
molecular bound-state energy Sb has been raised by 5 with respect to the atomic 
continuum. The difference in Zccman energy from the resonance position is also 
called the (experimental) detuning 6. As a result, the bare detuning 5b varies 
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Magnetic Field [G] Magnetic field [G] 

Figure 4: a) The scattering length as a function of magnetic-field strength near the very 
broad resonance of 6 Li at 834 G. It is given in terms of the Bohr radius ao- b) Hyperfine 
structure of the electronic ground state in 6 Li as a function of magnetic-field strength, where E 
is the internal energy and hp is Planck's constant. The nuclear spin of 6 Li is given by i a = 1, 
so that the spin of the outer electron gives rise to / a = 3/2 or / a = 1/2 for the total angular 
momentum. At large magnetic-field strength, the hyperfine energies are predominantly de- 
termined by the projection m B of the electron spin. The experiments with imbalanced spin 
mixtures that we consider later on in this review, are performed with states |1) and |2). 



linearly with the experimental detuning S, although they are not the same. As 
we may infer from Fig. [5] and as we also show next, they are separated by a 
constant shift. 

2.3. Dressed molecules 



We continue by looking for negative energy solutions of Eq. (Ill, meaning 



that we want to find the wavefunction and the energy of stable dressed molecules. 



We start by rewriting Eq. ( 11 ) to obtain the following equation 



(^m|K m G a (^)Km|^m> = E - S h , (12) 



with G a 1 (E) = E — Ho — Vb g . Assuming that we have been able to exactly solve 
the atomic part of the scatte 
we insert a completeness r( 
channel IV'k )> which gives 



the atomic part of the scattering problem, namely (Hq + Vbg)|'0k + ' 1 ) = 2ek|V'k + ' > )i 
we insert a completeness relation of the exact scattering states in the open 
/,(+)\ 



E - Sh = I j^r — e~2^ — = hm} > (13) 



where we have interpreted the expression containing the integral in Eq. ( 13 1 
as the molecular self-energy KE(E). Note that for simplicity we have ignored 
the possibility of bound states in the atomic background potential. If they are 
important, then they can be easily taken into account along the lines of Rcf. 
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Figure 5: Schematic representation of the physics near the resonance at Bo = 834 G of Li. 
The atomic continuum is in an electronic spin triplet channel, whose energy is lowered upon 
an increase of the magnetic field strength B. The bare molecule is in an electronic spin singlet 
channel and its bare energy is therefore not altered by the magnetic field. Due to interactions 
with the atomic continuum the bare molecule becomes dressed. When the dressed molecular 
energy, as sketched by the dashed-dotted line, reaches the atomic continuum the resonance 
takes place. In this section, we actually rotate the picture, namely the threshold of the atomic 
continuum stays at zero energy, while <5t, depends linearly on B. 



|120j . We thus see that the exact equation for the molecular bound state energy 
has the intuitive form E = 5b + KE(E). As we see soon, the resonance takes 
place when the energy E of the dressed molecule reaches the threshold of the 
atomic continuum, i.e. when E = and 5b — — 7iE(0). Therefore, we define 
the (experimental) detuning as 5 = 5b + fiS(O), so that we find for the energy 
equation E = 5 + KE'(E) with KS'(E) = KE(E) - 7iE(0). Since the resonance 
now takes place for 5 — 0, we indeed have that 5 = Afi(B — Bo). 



To make further progress with Eq. ( 13 ) we realize that the matrix element 



behaves for low momenta as 11171 



1 + ika hg ' ^ ^ 

where Ob g is the s-wave scattering length of the background interaction and g is 
by definition the coupling strength at zero momentum. In the limit of vanishing 
background interaction, the above equation reduces to (V'm|Kim|k) = <?, which 
means we have used a momentum-independent or local atom-molecule coupling. 
This approximation is appropriate, because the spatial extent of the molecular 
wavcfunction and the atom-molecule coupling is very small compared to the de 
Broglic wavelength of the scattering atoms. Using for the same reason also a 
local or momentum-independent background interaction Vbgj then the momen- 
tum dependence in Eq. |l4| ) is thus caused by the low-momentum behavior of 
the background transition matrix, which according to Eqs. ([7]), (JsJ) and ([9| is 
proportional to (1 + ifcabg) -1 . 
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Using Eqs. (131 and (14), we obtain that 



dk g 2 ( 1 1 



{2Tr) 3 l + k 2 al s \E-2e k 2e k 
1 — abgV —raE / H 



(15) 



with 77 = g 2 m 3 / 2 /4ttH 3 and where we have assumed a negative background scat- 
tering length, which is the case for 6 Li. The bound-state equation is now an 
analytically solvable cubic equation, although its solution is somewhat cumber- 
some [119] . A simpler result is obtained when we ignore a^ g , which is allowed 
close to resonance, where E goes to zero. Then, we find 



(16) 



so that for 5 — > we have E — — <5 2 /?/ 2 . This implies that the bound-state energy 
of the dressed molecule goes to zero quadratically with the applied magnetic- 
field strength B, although the bare molecular energy varies linearly with B. 




Moreover, from Eq. (Ill it follows that the dressed molecular state is given by 

IV^dr) = \fZ^\i>ra) + ^mlV'a), where 

= <V m |VUG a (£) 2 K m |V> m ) = -^fP. (17) 
z m oh, 



Eqs. (15) and (17) then allow us to calculate Z m (E), i.e. the amplitude of the 
dressed molecular state in the closed channel, where for the energy E we must 
substitute the solution of the bound-state equation. 

2. 4- Resonant atomic interaction 

Having determined the properties of the Feshbach molecules, we now ex- 
amine the effect of the Feshbach resonance on the atomic physics in the open 



channel. Solving Eq. (11) for the atomic channel, we hnd that 

{ffo + Vbg + Kn}!^) =£|Va), ( 18 ) 

where V m is the molecule-mediated interaction given by 

p _ Km|V'm)(V'm|Km (TQ\ 

Vm ~ E-S h ' [ ' 

and we also used that in our model |V'm)(V , m| is the unity matrix in the closed 
channel. As a result, the total interaction in the open channel is given by 
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V = Vb g + V m . From Eq. ([5|), we have that the exact transition matrix in the 
open channel is given by 



T = V + VGV 

EE f£+f£ (20) 



with G 1 (E) = E — Hq — Vbg — V m . The exact two-body T matrix is thus seen 

p2b 



to consist of a purely background part T?£ and a resonant part T 2 ^ . The first 



is given by T 2 ^ — Vb g + VbgG a Vb g , which, as before, leads for s-wave scattering 

-2b/ 



m. 



at zero momentum to T^(0) — AttK abg/ 

The momentum dependence of the resonant part may be further evaluated 

by 

(k |T res (2e k )|k) = ______ , (21) 



where we used Eq. (19 1. Moreover, we used that KE(E) is given by the left- 



hand side of Eq. (12), and that \ip m ) (ip m \ is the unity matrix in the molecular 



channel. Combining Eqs. ( |14[ ) and (21), we find that at zero momentum and 
energy T r 2 e ^(0) = -g 2 /(hY,(0) + 6 B ) = -g 2 /S. This indeed shows that for S = 
the transition matrix in the open channel goes to infinity at zero energy, implying 
a diverging scattering length. We thus have that 

T 2b( Q ) = 4 ^ 2 « _ 47r^ 2 a bg g 2 



m m 



6 



AB , 

!-^r . ( 22 ) 



fl bg 

m B - B 

where the width of the resonance AB, the location Bq and the background 
scattering length are readily determined experimentally by fitting to the data. 
As a result, the atom-molecule coupling constant g can be expressed in terms 
of the experimentally known parameters as g = v/ '^Trh^abgABAfi/m. 

To summarize, we have been able to solve for the properties of the dressed 
molecular state near a Feshbach resonance. Moreover, we have studied what the 
effect of this molecular state is on the transition matrix in the open channel. 
In particular, we have seen when the resonant part of the transition matrix 
diverges and gives rise to an infinite scattering length. Since the location of the 
resonance, its width, and the background scattering length are easily determined 
by fitting to the experimental data, we have access to all relevant parameters 



for a quantitative study of the molecular self-energy from Eq. ( 15 ) and the 



wavefunction rcnormalization factor Z m from Eq. (17). The results for Z m are 
shown in Fig. [6] for the broad resonance of 6 Li near Bq = 834 G, as measured 
in Ref. [TS] and as determined theoretically in Ref. [120] . It shows that the 
amplitude of the dressed molecule in the closed channel is very small over a 
very wide range near resonance. As a result, all the action happens in the open 
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Figure 6: Wavefunction renormalization factor Z m as a function of magnetic field strength, 
as adapted from Ref. 120 . Curve (a) is calculated without background interactions, curve 
(b) with a constant background scattering length a^g and (c) with a magnetic-field dependent 
ajjg as determined in Ref. [1211 . The experimental data (squares) is from Ref. |18| . At 
resonance, Z m goes to zero, while for positive detuning, the molecule is not stable in the 
two-body case. To get agreement with the experimental data also for positive detuning, a 
many-body calculation is needed |122| . 



channel, where, apart from allowing the resonance to actually take place, the 
bare molecular state hardly plays any role. This observation then leads to the 
so-called single-channel model of a Feshbach resonance, where only the atomic 
channel is taken into account with the transition matrix given by Eq. (22 1. The 
single-channel model is valid for wide resonances, where Z m is small. It is thus 
particularly valid for the extremely broad resonance of 6 Li near B§ = 834 G, 
which is used in the experiments on the imbalanced Fermi gas that are discussed 
later on. Therefore, we use the single-channel model in several theoretical many- 
body calculations throughout this review. Next, we use the single-channel model 
to qualitatively discuss an important application of Feshbach resonances, namely 
the study of the BEC-BCS crossover. 

2.5. Pairing in the BEC-BCS crossover 

If fcrmions of different spin attract each other, they have the tendency to 
form pairs, which may result in a paired condensate at ultralow temperatures 
[HIH]. However, if these fermions repel each other, then they have the tendency 
to align in spin space, favoring (itinerant) ferromagnetism [123, 1241 1125j . Al- 
though close to resonance these two instabilities compete [126] . we consider in 
this review solely pair condensation. Since in ultracold atomic gases the range 
of the interaction is typically very small compared to the interparticle distance, 
it is convenient to use a contact potential of strength Vg as interaction poten- 
tial V(r), i.e., V(r) — Vo<5(r), to incorporate interaction effects. The Fourier 
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Figure 7: Schematic representation of the BEC-BCS crossover. In the BCS regime, the 
microscopic interaction strength, given by Vo, and the effective interacting strength, given 
by the scattering length a, are both negative and small, namely |o| < n -1 / 3 with n the 
total particle density. Pairing is then a many-body effect and the size of the Cooper pairs is 
much larger than the average interparticle distance. When \a\ > n~ 1 / 3 , we enter the strongly 
interacting regime. At Vo = — 2n 2 h 2 /mAo (see text), the scattering length diverges, which is 
also called the unitarity limit. Here, the pair size is comparable to the average interparticle 
distance. For stronger microscopic attractions, a two-body molecular bound state enters the 
interaction potential. When \a\ < n — we enter the BEC regime. The size of the molecules 
is here much smaller than the average interparticle distance and the ground state of the system 
is a weakly interacting molecular Bose gas. 



transform of the potential is then a constant, which leads to a particularly sim- 
ple form of the Lippmann-Schwinger equation in momentum space. For s-wave 
collisions, the transition matrix at zero momentum and energy becomes [30] 



1 1 



T 2b (0) Vo V ^ 2e k 



where the kinetic energy of a particle with mass m and wavevector k equals 
ek = h 2 k 2 /2m and V is the volume of the systemj^] Note that the sum on the 



right-hand side of Eq. (23) is actually not convergent, which comes from the 
anomalous behavior of the contact potential at high momenta. This is usually 
not a problem, since we are only interested in the low-energy properties of 



the interaction. As we see explicitly in Section 4.2 we can then use Eq. (23 1 
to eliminate the contact potential from our theory in favor of the two-body 
T matrix, which at zero energy is related to the s-wave scattering length a 
through T 2b (0) = 4Trh 2 a/m. This scattering length can be directly extracted 
from experiment. 

Another approach to regularize the high-momentum behavior of the contact 



1 ln actual analytical or numerical evaluations of the sum over states in k-space, we make 
use of the continuum limit f/V — > J dk/(2-7r) 3 . 
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potential is to consider only momenta up to a certain cut-off Ao for the sum in 
Eq. ((23]), after which we find [T27] 

4irh 2 a 7r . 
Vo = ^-r-, 24 

whose result is plotted in Fig. [7] The weakly interacting regime, where both a 
and Vo are small and negative, is also called the BCS regime, since it is analogous 
to the weakly attractive case known from ordinary superconductors. In the BCS 
regime, the two-body potential is too weak to support a bound state, so that 
Cooper pairing is truly a many-body effect that occurs only in the presence of a 
Fermi sea [7] . The size of these Cooper pairs then turns out to be much larger 
than the average interparticle distance. When Vq = —2TT 2 h 2 /mA , we see that 
the scattering length diverges. This resonance in the cross section is physically 
caused by a bound state that enters the attractive potential |128j . To see the 
latter more clearly, we can consider a slightly more realistic interaction potential, 
namely a square well. The exact solution to the resulting problem shows that, 
upon increase of the well depth, the scattering length diverges each time a new 
bound state enters the square- well potential |117j . Such resonances are called 
single-channel shape resonances. Note that it is not so clear how to tune a shape 
resonance experimentally, because it is usually not possible to precisely control 
the interatomic interaction potential. However, in the previous paragraphs, we 
discussed the more flexible mechanism of a two-channel Feshbach resonance, 
which is nowadays routinely used by experimentalists to vary the scattering 
length at will. 

The region where the scattering length becomes infinite is commonly referred 
to as the unitarity limit. Here, the size of the Cooper pairs turns out to be 
comparable to the average interparticle distance. The unitary regime is also 
called the strongly interacting limit. This as opposed to the regime where the 
microscopic attraction has become so strong that there is a deep bound state 
in the interaction potential. As a result, bosonic molecules are formed, whose 
size is much smaller than the average interparticle distance. This regime is 
then called the BEC regime. It is unique that the complete evolution from the 
BEC regime to the BCS regime can be explored in ultracold quantum gases 
[ID]. From Fig. [7] we might be surprised that the evolution merely leads to a 
crossover. Namely, the system is seen to evolve through a true resonance in 
the effective interaction strength, which we might have expected to profoundly 
influence the thermodynamics of the mixture. Moreover, theoretical predictions 
based on perturbative approaches are expected to fail near resonance. This 
particularly holds for the mean-field BCS theory, which was used to predict the 
crossover [13] . As a result, there has been much theoretical and experimental 
interest in the BEC-BCS crossover, and in particular on the behavior of the 
Fermi mixture in the unitary regime [101 129] . 

A second reason for the appeal of the crossover is the unification of BEC- 
like superfluidity and BCS-like superfluidity as two-sides of the same coin. Both 
can be viewed as a coherent state (or Bose-Einstein condensate) of fermionic 
pairs, which is signalled by a nonzero expectation value of the pair annihilation 
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operator, i.e. (?/>_(R + r/2)?/>+(R— r/2)). Here, we introduced the annihilation 
operator ?/v( x ) for a single atom with spin a at position x, while R is the 
center-of-mass coordinate and r the relative coordinate of the pair. This puts 
us in the position to define the equilibrium BCS order parameter (A), given by 

(A)(R) = |rfrF(r)(v;_(R+|)^ + (R-|)) 
= Vb^-(RhMR)) 

= 7E(i J f + rf + ,!-^' , ' R (25) 

q,k 

where fiq is the center-of-mass momentum and Kk is the relative momentum. 
Because the equation above is seen to involve an integral over the attractive 
interaction, the order parameter also describes the energy cost to break up a 
Cooper pair. For this reason, (A)(R) is also referred to as the (local) pairing- 
gap. Usually, the paired state only occurs at at zero momentum, hq = 0, 
which upon pair condensation gives rise to an order parameter that is spatially 
independent, namely 

< A > = (26) 

k 

The above equation also reveals that in this case the Cooper pairing occurs in 
momentum space between particles of opposite spin and momentum. 

Using BCS mean-field theory it is only possible to study the crossover qual- 
itatively at zero temperature, where it describes the evolution from a conden- 
sate of loosely-bound Cooper pairs to a condensate of tightly bound molecules 
OH!]. For the critical temperature curve, BCS mean-field theory fails even on 
a qualitative level. This is because the BCS critical temperature describes physi- 
cally the effect of pair-breaking, which is not the correct mechanism for the phase 
transition in the BEC regime, where superfluidity is lost to the thermal occupa- 
tion of nonzero momentum states by tightly-bound pairs. The simplest theory 
that also describes these pairs with nonzero momentum is called the Nozieres- 
Schmitt-Rink approximation 129 , which we discuss in more detail in Section 
|4.2| It turns out to give a rather accurate description of the critical temperature 
throughout the crossover [130J. Although the BEC-BCS crossover was predicted 
already a long time ago, it was hard to experimentally verify the smooth evolu- 
tion from BCS to BEC superfluidity using condensed-matter systems, since the 
attractive phonon-mediated electron-electron interaction is not so easily tun- 
able. Therefore, as mentioned in the Introduction, the BEC-BCS crossover has 
only been observed experimentally in recent years [MJ El [TBI H3 02] , profiting 
from the extreme tunability of ultracold quantum gases. The experiments on 
imbalanced Fermi gases that we discuss more extensively in this review were 
performed mainly about halfway the crossover, where the scattering length di- 
verges. This is also called the unitarity regime, where the Fermi gas has an 
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exceptionally large critical temperature that is about one-tenth of the Fermi en- 
ergy [TH], and where the Fermi gas also shows remarkable universal properties 
pTl[T82irT33llT34] . 

3. Mean-field theory 

In this section, we use BCS mean-field theory to study the strongly inter- 
acting two-component Fermi mixture with a population imbalance. We start 
out by giving a comprehensive historical overview of the experiments that have 
been performed since Zwierlein et al. [24] and Partridge et al. [25] achieved 
control over the spin polarization of the Fermi gas. After this experimental 
overview, we calculate the mean-field phase diagram for the homogeneous case. 
Although pcrturbative approaches are not expected to be quantitatively correct 
for strong interactions, mean-field theory turns out to be useful for explaining 
the relevant physics in the system. In order to also arrive at the phase diagram 
for the trapped case, we use the local-density approximation, which assumes 
that the gas behaves locally in the trap as if it were homogeneous. From the 
mean-field phase diagram, we can understand the qualitative aspects of the 
early experiments exploring the imbalanced Fermi mixture. We end this section 
on mean-field theory with a calculation of the phase diagram of the 6 Li- 40 K 
mixture. 

3.1. Experimental overview 

In the beginning of 2006, two experimental groups obtained a full control 
over the spin imbalance or polarization P in an ultracold atomic Fermi mixture. 
This polarization is defined by P = (N + - N-)/(N + + AT_) with N a the total 
number of particles in the hyperfine state a. These first experiments at MIT 
by Zwierlein et al. [21] and at Rice University by Partridge et al. [35] realized a 
strongly interacting mixture of 6 Li atoms in the two lowest hyperfine states of 
the ground state (see Fig. [4]), where the atom number in each of the two states 
could be precisely tuned. Within a few months time these experiments were 
followed by a large amount of theoretical studies, which studied the imbalanced 
Fermi gas using mean-field approaches |1001 11351 11361 11371 11381 |8"S] , using ther- 
modynamic approaches |139| I140j , using diagrammatic approaches |141) , using 
the Bogoliubov-de Gennes approach }142[ 1143] , including the effects of surface 
tension [144] . including the gradient energy of the density profiles |145] . and 
including Gaussian order parameter fluctuations |146| . The interesting aspect 
of a polarization in the system is that not all atoms can find a partner to pair 
up with, because the fermions only have attractive s-wave interactions when 
they are in different hyperfine states. Since pairing is the mechanism behind 
fermionic superfluidity, the following fundamental question arises: what hap- 
pens to a superfluid paired Fermi mixture upon increasing the polarization? 

The early 6 Li experiments of Zwierlein et al. [Ml H47j and Partridge et al. 
|25| were performed in a trap to confine the atomic clouds in space. Both exper- 
iments revealed that superfluidity in an ultracold Fermi mixture with attractive 
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Figure 8: Sketchy impression of three different phases that can be present for a strongly 
interacting Fermi mixture in the trap. The x and y coordinates correspond to the axial and 
the radial direction in the trap, while the z coordinate represents the total density of atoms. 
Although the axial direction is elongated in the actual experiments, we have drawn here a 
more spherically symmetric situation, a) At high temperatures, the gas is in the normal 
phase (N), leading to a thermal distribution of atoms in the trap. b)-c) At low temperatures, 
a condensate of Cooper pairs appears in the center of the trap, where the particle densities 
are highest. The condensate leads to a pronounced enhancement of the central densities. The 
superfluid core and the normal outer region in the trap can be separated by b) a continuous 
second-order transition (S) or c) a discontinuous first-order transition (P). Note that we have 
exaggerated the 'bump' in the profile due to the Cooper-pair condensate and the jump in the 
density due to the first-order phase transition for illustrative reasons. 



interactions is in first instance maintained upon going to an imbalance in spin 
populations. However, on a more detailed level contradictory results were ob- 
tained for the behavior of the trapped mixture as a function of polarization. 
Namely, Zwierlein et al. observed a rather smooth phase transition between a 
phase with a superfluid core and a phase that was fully normal at a high crit- 
ical polarization of about 0.7. However, Partridge et al. seemed to observe a 
transition between two different trapped superfluid phases at a low critical po- 
larization of about 0.1. Moreover, Partridge et al. did not observe a vanishing of 
the superfluid core in their experiments, where even at their highest imbalances 
(P > 0.90) the core seemed fully paired. 

More experimental and theoretical studies followed in order to try to un- 
derstand these differences. In a later study at MIT, Shin et al. showed that 
their trapped Fermi mixture was described by a shell structure consisting of an 
equal-density superfluid core, surrounded by a partially polarized normal shell 
and a fully polarized outer region [148j . As a function of polarization, the MIT 
data showed a smooth onset of the superfluid condensate in the center of the 
trap, which agreed with a continuous, or second-order, phase transition. Also at 
Rice University a further study was performed, where Partridge et al. observed 
a deformation of the superfluid core at their lowest temperatures [71 . This 
deformation was explained in terms of a first-order interface between the super- 
fluid core and the normal outer region in the trap [1441 1149] . see also Fig. [8jc). 
Namely, a first-order interface in general gives rise to a surface tension, which 
energetically favors a minimal area of the interface, i.e. a spherical shape. How- 



ever, the trap that was used in the experiments at Rice university was actually 
highly elongated, so that the superfluid core was consequently deformed from 
the trap shape by the surface tension. At higher temperatures the deformation 
disappeared, although the core seemed to remain paired. The natural expla- 
nation of this behavior is that at higher temperatures the superfluid-normal 
transition in the trap is of second order, as in Fig. [sjb) , so that there is no 
surface tension and deformation. 

Since the higher-temperature results at Rice University resembled qualita- 
tively the behavior that was seen in the MIT experiments, it seemed like a 
difference in temperature would be the most natural way to explain both exper- 
iments in a single theoretical picture |85j . Moreover, in another experiment at 
MIT, Shin et al. performed precise measurements as a function of position in the 
trap and temperature [72] ■ Because the local-density approximation applies for 
the MIT experiments, this actually means that Shin et al. were able to experi- 
mentally map out the homogeneous phase diagram of the spin-imbalanced Fermi 
mixture in the unitarity limit. They now observed both a second-order phase 
transition in the trap at higher temperatures as well as a first-order transition at 
their lowest temperatures. On a qualitative level, all experimental results then 
seem to fit in a phase diagram that is governed by a tricritical point [85] [88] . 

However, on the quantitative level there remained striking differences be- 
tween the two experimental groups. First of all, there was the difference in 
the critical imbalance P c , at which the trapped gas becomes completely nor- 
mal. Namely, Shin et al. found at their lowest temperatures P c < 0.8 [1471 172"] . 
while Partridge et al. found P c > 0.9 [71]. Although the first result agrees 
with Monte-Carlo calculations in the local-density approximation (LDA) [150] . 
the latter result was not understood. The second issue involves the normal 
region that surrounds the superfluid core. While the MIT experiments found 
that the normal state is partially polarized close to the superfluid interface, the 
Rice experiments only observed a fully polarized normal state. The absence 
of the partially polarized normal region in the Rice experiment was also not 
understood. The last difference is the observation of deformation by the Rice 
experiment, whereas in the MIT experiment such deformations were not seen. 
The high value of the surface tension that is needed to explain the deformation 
observed at Rice University could also not be microscopically accounted for. 

More recently, the imbalanced Fermi gas was also experimentally realized 
at ENS by Nascimbene et al., who found very similar results to the MIT ex- 
periments [73[ O E5] . The latest experiments by the group at Rice indicate 
that their deviating experimental results are caused by spin currents in the 
trap that create a long lived metastable state with a deformed superfluid core 
|78j . These experiments confirm the theoretical proposal of Parish and Husc 
|151j . who showed that this metastable state could arise due to the elongated 
trap geometry, in which evaporation occurs predominantly from the trap center. 
Once phase separation occurs in this geometry, particle and heat transport are 
suppressed at the narrow interface between the gapped superfluid core and the 
normal outer region |1521 11531 1151] . It turns out that the resulting evapora- 
tive depolarization of the central core stabilizes superfluidity, and thus favors 
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the superfluid phase over the normal phase, even when the overall polariza- 
tion of the gas is above the expected critical imbalance. The suppression of 
transport across the interface allows the non-equilibrium density distribution to 
remain metastable with observed lifetimes of multiple seconds [78]. Although 
spin transport is an interesting topic that is currently actively being studied 
experimentally for imbalanced Fermi gases [751 177] , we treat in this review only 
the equilibrium case. This also means that in the following we mainly com- 
pare the theoretical results with the experiments performed initially at MIT 
|231 H571 11551 172]. 

3.2. Mean- field grand potential 

We start our theoretical discussion of the strongly interacting Fermi mixture 
at the mean-field level, which is useful for qualitatively explaining the relevant 
physics. The derivation of the grand-canonical thermodynamic potential density 
that follows from mean-field BCS theory can be found in several textbooks, see 
for example Ref. |154j using operator methods or Ref. [30] using functional 
methods. In Ref. [89] , the case for both a spin- and mass- imbalance is explicitly 
derived. It yields for the solely population-imbalanced case at unitarity that 

1 f |A| 2 1 
^bcs [A; T,n„\ = y ^ | e k — ^ — ^o; k + | 

-^EMl + e-^), (27) 

cr,k 

where the kinetic energy of the atoms is given by ek = h 2 ^ 2 /2m, m is the mass 
of the fermions, V is the volume, T is the temperature and f3 — 1/k^T. The 
physical meaning of the BCS order parameter A was explained in section |2.5| 
The index a = ± specifies the hyperfine state and is also called the (pseudo)spin 



of the fermions. There are two slight differences between Eq. ( 27 1 and the stan- 
dard BCS expression for the grand-canonical thermodynamic potential. First of 
all, there is a term — |A| 2 /T 2b missing, where the two-body transition matrix is 
given by T 2b = 4irh 2 a/m with a the s-wave scattering length. This is because 
the scattering length diverges in the unitarity limit, which is the regime where 
the experiments operate. Second, we allow the chemical potentials for the two 
hyperfine states [i& to be unequal, since this takes into account the population 
imbalance in the two spin species that is realized by the experiments. The aver- 
age chemical potential fi is given by fi = + /i_)/2, while half the difference 
is denoted as h = — /u_)/2. In the case of nonzero h, the dispersions of 
the Bogoliubov quasiparticles fto; CTi k are spin-dependent. Namely, we have that 
^w CT ,k = —ah + fi^k with hujk = yj (ek — A*) 2 + I A| 2 |155) . This follows from 
the usual Bogoliubov diagonalization of the mean- field hamiltonian |156j . The 



logarithms in Eq. (27) describe an ideal gas of fermionic quasiparticles with 
dispersion fiw^k, while the other terms describe the equal-density Cooper pair 
condensate. 
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Figure 9: a) The behavior of the grand-canonical thermodynamic potential density wbcs 
as a function of the BCS order parameter [A| for a second-order phase transition. For high 
temperatures T, the global minimum is at |(A}| = 0, i.e. the normal state. At the critical 
temperature T c the minimum becomes a maximum. For T < T c the system is in the superfluid 
state, b) Corresponding behavior of the expectation value of the order parameter |{A)|. At 
T c , | (A) | goes continuously to zero, c) Same as a), but now for a first-order transition. At the 
critical temperature T c , the separated minima of the normal and the superfluid state have the 
same pressure, given by p g = — ojbcs> an d |(A}| makes a jump. This discontinuity is shown 
in panel d). 



To determine the atomic density n a in spin state a, we use the relation 
n a = — <9wbcs[A; T, /vj/ctyv, which leads to 

n CT [A;T,/v] = 1 {K| 2 /(?^,k) + K| 2 [l - /(^- CT ,k)]} , (28) 

k 

where f(x) = l/(e l3x +l) are the Fermi distributions. The BCS coherence factors 
|Mk| 2 and |^k| 2 are determined by the relations |wk| 2 = [1 + (ek — A*)/^k]/2 and 

Ki 2 + kki 2 = i mm- 

Depending on the temperature T and the chemical potentials the ther- 
modynamic potential density wbcs can give rise to either one, two or three 
extremal points. At high temperatures there is only one extremum, namely a 
global minimum at A = 0, so that the system is in the normal phase. For the 
well-studied balanced case, h = 0, BCS theory predicts that at a certain critical 
temperature T c the extremum at A = becomes a local maximum. The global 
minimum then continuously shifts away to a nonzero value of A, and the system 
enters the superfluid phase. Since the transition evolves continuously as a func- 
tion of temperature, it is called a continuous, or second-order, phase transition 
and it is shown in Figs, una) and (b). For h ^ 0, the minimum at A = can 
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a) 



b) 



Figure 10: Quasiparticle dispersions for a) the BCS superfluid phase and b) the Sarma 
superfluid phase. The upper branch gives the dispersion for the spin-down quasiparticles, 
that consist of spin-down particles and spin-up holes. The lower branch gives the dispersion 
for the spin-up quasiparticles, that consist of spin-up particles and spin-down holes. These 
dispersions have their minima given by |A| ± h at wavevectors for which = fi. In the 
BCS case the quasiparticle spectra are gapped. In the Sarma case, a part of the spin-up 
quasiparticle branch is below zero, such that its filling lowers the ground-state energy. As a 
result, additional spin-down holes and spin-up particles enter the ground state leading to a 
polarized superfluid. 



also be a local minimum, so that there is both a local maximum and a global 
minimum at values of A unequal to zero. As is seen in Figs. [9^c) and (d), this 
can cause a discontinuous, or first-order, phase transition. The extrema of the 
thermodynamic potential density can be found by differentiating with respect 
to A* and equating the result to zero. As the above discussion implies, there 
is always one solution given by A = 0. The other solutions are found from the 
so-called BCS gap equation 



v 4^ 







(29) 



which thus has either one or two solutions. The study of the extrema of the 
thermodynamic potential allows for a determination of the phase diagram as a 
function of the chemical potentials and the temperature, which we perform in 
Section l3~4l 



3.3. Sarma phase 

But first, let us briefly discuss in more detail the homogeneous superfluid 
phases that we encounter in the spin-imbalanced case. Below the critical tem- 
perature T c , we have that |A| ^ 0, in which case we distinguish between two 
possibilities. Namely, we have either that h < |A|, or that h > |A|. The first 
case we call a BCS superfluid, because, as we see next, it corresponds to the 
fully-gapped situation known from ordinary BCS superconductivity in metals 
[2 . The second case leads to a so-called Sarma superfluid, which gives rise to 
a gapless quasiparticle dispersion for the majority spin species fiwk,+ , as was 
first discussed by Sarma [84]. The Sarma phase is sometimes also referred to 
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as the interior-gap phase or the breached-pair phase |55J [EE] ■ To see the differ- 
ence between the two phases more clearly, we discuss the notion of a Bogoliubov 
quasiparticle and the meaning of its dispersion in more detail. A spin-up (down) 
quasi-particle is a linear combination of a spin-up (down) particle and a spin- 
down (up) hole, see e.g. |157j . namely 

&+ )k = u ^+m ~ "kfl- ,-k, (30) 
b\ _ = Wkfl+,k + v k a[_ _ k , (31) 

where k creates a particle with wavevector k and spin a, while 6^ k creates a 
quasiparticle. Note that we have chosen real coherence factors Mk and i>k- Quasi- 
particles behave like ordinary particles in the sense that they can be assigned 
an energy, wavevector and spin. 

Physically, a quasiparticle excitation describes a single-particle excitation on 
top of the Cooper pair condensate. To make such an excitation, first a Cooper 
pair has to be separated into two uncorrelated atoms, after which one of the 
particles is taken out of the system, while the other remains. The remaining 
particle then determines the spin of the quasiparticle excitation. The energy 
difference between the excitations of different spin is thus 2h, namely the differ- 
ence in the chemical potentials between the two spin states. The quasiparticle 



dispersions are shown in Fig. 10 where panel (a) corresponds to the BCS case 



and panel (b) to the Sarma case. The upper curve in Fig. 10 'a) shows the 
spin-down quasiparticle spectrum huj-^ and the lower curve shows the spin- up 
quasiparticle spectrum huj + ^. It is seen that it cost a certain nonzero amount 
of energy to make an excitation in one of the two branches, which are therefore 
said to be fully gapped. At zero temperature, both branches are completely 
empty, so that only the equal-density BCS ground state, describing the conden- 
sate of Cooper pairs, remains. In Fig. |10| h), the same quasiparticle dispersions 
are drawn for the gapless Sarma phase, which arises when h > |A|. The curves 
have the same interpretation as for the gapped BCS case, however, now the 
branch for the majority spin-up quasiparticles fiw+,k goes through zero at the 
wavevectors k satisfying £k = M± \/h 2 — |A| 2 . As a result, this additional part 
of the quasiparticle branch below zero is filled at zero temperature, because in 
this way the ground-state energy is lowered. The spin-up quasiparticles make 
the Cooper pairs vanish in the corresponding momentum range, so that here 
only spin-up particles remain and the superfluid becomes polarized. 

We see this more clearly if we explicitly calculate the average occupation 
number for a single-particle quantum state that is specified by the wavevector 
k and spin a. We designate this occupation number with N a (ls.) and it is given 



by the expression that is summed over in Eq. (28), namely 

N a (k) = |u k | 2 /(^ CT ,k) + M 2 [l - /(fcj-„,k)]. (32) 

Indeed, summation over all occupation numbers A CT (k) yields the total number 
of particles in state a. Due to the Pauli principle, A CT (k) can be maximally 



unity. In Fig. 11 we show the result for N a (k) both in the case of a BCS su- 



perfluid and a Sarma superfluid at zero temperature. In the BCS case of Fig. 
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Figure 11: Average occupation numbers 7V CT (k) of single-particle quantum states with mo- 
mentum k and spin a for the Sarma phase. At wavevectors for which < fi — ^Jh 2 — |A| 2 
or e k > A 1 + \/h 2 — |A| 2 , A^k) equals v 2 . corresponding to the occupation numbers of an 
equal-density BCS superfluid. For wavevectors in between, ./V_|_(k) = 1 (dashed line) and 
iV_(k) = (dash-dotted line), so that the superfluid becomes polarized. 



[lO[a), we obtain from Eq. (28 1 that the occupation for both spin states is given 
by -/V±(k) = v£ & t zero temperature, so that the particle numbers are the same 



and the superfluid is fully balanced. However, in the Sarma case of Fig. 10 b) , 
the additional spin-up quasiparticles, which create extra spin-down holes and 
spin-up particles, are seen to cause a full polarization near the average chemical 
potential /i. For the wavevectors leading to a positive quasiparticle dispersion 
Sw +i k, the BCS behavior is recovered. The resulting configuration for the Sarma 
phase is sometimes also referred to as 'phase separation' in momentum space, 
because for small |k| and for large |k| the system has the supcrfluid-state occu- 
pation numbers from BCS theory, while for |k| around {2ra^i) 1 / 2 /h the states are 
occupied as a fully polarized normal state. Fig. [TT] shows that the Sarma phase 
can also be defined at zero temperature as a superfluid with a Fermi surface, or 
with multiple Fermi surfaces. At nonzero temperatures, the sudden rising and 



lowering of the occupation numbers in Fig. 11 becomes smoother |158j . until at 
high temperatures this nonmonotonic behavior completely disappears. 

If the Sarma phase were stable at zero temperature, it would be separated 
from the BCS phase by a true quantum phase transition, meaning that non- 
analytic behavior in thermodynamic quantities would be observed as the ma- 
jority quasiparticle spectrum became gapless. However, as we find in the next 
section, this gapped-gapless transition is preempted by a first-order phase tran- 
sition to the normal state, so that the situation |A| > h is found to be never 
stable at unitarity. More precisely, if we study at zero tempertaure the behavior 
of the equilibrium order parameter (A), corresponding to the global minimum 
of the thermodynamic potential, then we find that when |(A)| = 0.70/i, a first- 
order transition to the normal state takes place, where the latter remains the 
global minimum for larger h. However, at nonzero temperatures, the situation 
is very different. Namely, close to a second-order phase transition (A) becomes 
arbitrarily small, and the condition for stable gapless superfluidity |(A)| < h is 
readily satisfied. Moreover, in the case of a large mass imbalance between the 
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fermionic species, the Breached-Pair 1 phase (Sarma phase with a single Fermi 
surface) is expected to be stable all the way down to zero temperature [85] . The 
same holds for the BEC side of the resonance where the Sarma or Breached-Pair 
1 phase is stable at zero temperature according to mean- field theory |100| . The 
sharp Fermi surfaces in the momentum distribution of Fig. [TT] can be used to 
detect the Sarma phase experimentally and distinguish it from the BCS phase. 
A precise scheme to measure such momentum distributions through time-of 
flight Raman imaging was proposed by Yi and Duan |158j . At nonzero tem- 
peratures, the Fermi surfaces of the Sarma phase are smoothened, so that the 
evolution from the gapped to the gapless regime is no longer a phase transition, 
but rather turns into a smooth crossover. As a result, there is no non-analytic 
thermodynamic behavior during this evolution anymore. The sharp Fermi sur- 
faces have now become continuous, but still there are characteristic bumps and 
dips in the momentum distributions. These non-monotonic distributions could 
be a clear experimental sign of being in the Sarma regime [158] . Whether or 
not this non-monotonic behavior can be observed experimentally at unitarity 
for the mass-balanced case is an open question, both theoretically and experi- 
mentally. For the mass-imbalanced case and on the BEC side of the resonance, 
the Sarma phase is stable on the mean-field level down to very low, even zero, 
temperature, resulting in observable Fermi surfaces. 

3.4- Homogeneous phase diagram 

Having looked in more detail at two homogeneous superfluid phases that can 
occur in the imbalanced system, we are now in the position to determine the 
phase diagram for the two-component Fermi mixture in the unitarity limit. In 
first instance, we calculate this diagram as a function of the temperature and 
the chemical potentials. This means that for any given combination of T and 
fi a , we specify whether the global minimum of wbcs is realized at a zero or a 
nonzero order parameter (A) . This determines whether we are in a superfluid or 
in a normal phase. For the superfluid phase we also specify whether h is larger 
than | (A) | or not. The first case then leads to the gapless Sarma regime, while 
the second leads to the gapped BCS regime. To calculate the phase diagram we 
thus need conditions that specify the phase boundary between the superfluid 
and normal phases, for which there are two possibilities. Namely the transition 
can be either continuous or discontinuous. This was illustrated in Fig. [9] 

For the continuous phase transition, the condition is given by the minimum 
at A = becoming a maximum, or equivalently, the second derivative of wbcs 
changing sign from positive to negative. The criterion thus becomes 
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where we introduced the shorthand notation /<x(k) = l/[e@( ek ^ + 1]. Note 
that this condition is only valid when the minimum of cjbcs at | A| = is a global 
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minimum before turning into a maximum, as in Fig. [9|a). This is something 
that we have to check. It is certainly the case if, at the critical temperature T c , 
cjbcs only has nonnegative and positive coefficients for its expansion in |A| 2 , 
which is true for the balanced case, h = 0. Following the line of second-order 



phase transitions for increasing h/fi, see the solid line in Fig. 12 a), then the 
first higher-order coefficient of wbcs that turns negative is the fourth-order one, 
which is proportional to 
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When both cth — = 0, then we are at the so-called tricritical point (TCP). 
Here, the continuous phase transition turns into a discontinuous one. Below 
the tricritical point, the condition for the first-order transition line is given by 
the disconnected normal and superfluid minima having an equal grand-potential 
density, or cj B cs[0; T c , fi a ] = w B cs[(A);T c , /v]- This is illustrated in Fig. ^c), 
where the corresponding behavior of the order parameter as a function of tem- 
perature is sketched in panel (d). 

The conditions for the phase boundaries give rise to the phase diagram of 



Fig. 12 a). Also drawn is the crossover between gapped BCS and gapless Sarma 
superfluidity, given by the condition |(A)| = h. A special feature of the phase 
diagram is that it has a certain kind of universality, due to the divergence of 
the scattering length |131j . Namely, because now the interaction strength does 
not provide a physical energy scale, the only energy scales left are associated 
with the temperature and the particle densities in the system. The latter scales 
can be characterized by the chemical potentials or the Fermi energies epo- 



h 2 (6n 2 n a ) 2 / 3 /2m. By scaling all energies in Fig. 12 'a) with the average chemical 
potential fj,, the phase diagram becomes only dependent on T/fi and h/ \x and 
not on the specific value of fj, anymore. At zero temperature, the first-order 
superfluid-normal transition occurs when h c = 0.81/x, at which the gap jumps 
from | (A) | = 1.16/i to zero. This critical value for h is sometimes also called the 
Chandrasekhar-Clogston limit [159L 1160] . We thus have that |(A)| > h c , so that 
the transition is between an equal density superfluid and a polarized normal 
state. This means that according to mean-field theory, the Sarma phase does 
not occur at zero temperature. 

In the normal region of the phase diagram, we could in principle make a 
further distinction between two different normal phases, namely a partially po- 
larized normal phase, where the minority species is still present, and a fully 
polarized normal state, where only the majority species remains. At zero tem- 
perature, the transition between this partially polarized phase and the fully 
polarized phase happens when h = /i in mean-field theory. In the phase dia- 
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Figure 12: a) The phase diagram of the homogeneous two-component Fermi mixture in the 
unitarity limit, consisting of the gaplcss supcrfluid Sarma phase (S), the gapped supcrfluid 
BCS phase and the normal phase (N). The transition from supcrfluid to normal can be cither 
continuous (full line) or discontinuous (dashed line), and the two possibilities meet at the 
tricritical point (TCP). Between the BCS regime and the Sarma regime of superfluidity there 
is a crossover (dash-dotted line). Both the temperature T and half the chemical potential 
difference h are scaled with the average chemical potential fi. b) The same diagram but now 
as a function of the polarization p = (n+ — n_)/(n+ + n_) and with the temperature scaled 
by the Fermi temperature of the majority species Tp + = ep + /kg. Due to the discontinuous 
nature of the transition below the tricritical point there is a jump in the polarization, causing a 
forbidden region (FR) in the phase diagram where the gas is unstable against phase separation. 
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grams that follow we typically do not explicitly make the distinction between 
the two cases, so that we refer to both phases together simply as the normal 
phase. 

The first-order transition line in Fig. [l2^a) is characterized by a jump in the 



order parameter (A). From Eq. (28), we find that this causes also a jump in 
the particle densities and in the polarization p = (n + — n_ ) / (n + + n_ ) . At zero 
temperature the discontinuity in polarization is largest and mean-field theory 
predicts a jump from an unpolarized superfluid, p — 0, to a normal state with 
a critical polarization of p c = 0.93. When we calculate for each point in the 
phase diagram of Fig. [l2^a) the corresponding particle densities, we find the 
diagram of Fig. 12 'b). The phase boundaries are given as a function of p and 
T/Tp + , where the Fermi temperature is defined by fceTjv = ep a . The phase 
diagram is again universal, so that the phases are uniquely determined by the 
polarization p and the ratio T/Tp a . The diagram in Fig. [l2|b) is seen to have 
a forbidden region (FR), which means that the combinations of temperature 
and polarization inside that region are not stable. At zero temperature, we see 
for example that the whole region between p = and p = 0.93 is forbidden. 
However, as mentioned in the introduction to this chapter, experimentally the 
polarization can be fixed at any value. If the system is forced to be in the forbid- 
den region at a certain polarization pf, then phase separation occurs (21)1 1161] . 
It means that the system forms both superfluid domains with low polarization 
and normal domains with high polarization, so that in total p = pi is satisfied. 
We end this discussion by noting that the present mean-field study so far only 
considers homogeneous s-wave superfluid phases. There have been theoretical 
proposals claiming that small parts of the phase diagram at unitarity are occu- 
pied by more exotic superfluid phases, such as FF or LO-like superfluids [M] . 
or induced p-wave superfluids These more exotic possibilities are not in- 
cluded by the present mean-field study. So far, they have also not yet been seen 
in experiments. 



3.5. Local-density approximation 

All experimental set-ups for ultracold atomic quantum gases invoke an ex- 
ternal trapping potential, which is needed to keep the atom cloud together for 
the duration of the experiment. The trapping potential is created by an ex- 
ternal magnetic field, that acts on the magnetic dipole moment of the atoms, 
or by the strong electric field in a laser beam, which induces an electric dipole 
moment in the atoms. As a result, inhomogeneous magnetic or electric fields 
become potential energy landscapes for the atoms, with which the particles can 
be confined in space without using material walls. Since the cold atoms accu- 
mulate around the minimum of the potential energy, the trap can typically be 
well approximated by a harmonic potential. Thus, if we want to describe an 
actual quantum gas experiment, we inevitably must study the effect of the ex- 
ternal potential. This can be somewhat inconvenient, because inhomogeneous 
systems are typically more cumbersome to deal with theoretically. Fortunately, 
if the trapping potential is locally flat enough, then we may consider the gas 
to be homogeneous at that point in the trap. This is the physical essence of 
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Figure 13: The homogeneous phase diagram of Fig. |12| but now with 6 arrows (a)-(f) that 
represent 6 different polarized phases in a trap. The tail of the arrow represents the center of 
the trap, while the decrease of fi(r) for increasing radius r causes the arrow to move up and 
to the right, leading the arrow to the normal phase. For each arrow, we also draw a schematic 
representation of the corresponding trapped phase. The center of the spheres correspond to 
the centers of the trap, while the thick full line denotes the edge of the cloud. The core of the 
cloud can be either in the normal phase (N), the gapped superfluid phase (BCS) or the gapless 
superfluid phase (S). These locally homogeneous phases are separated in the trap by either 
a second-order transition (thin full line), a first-order transition (dashed line) or a crossover 
(dash-dotted line). 



the local-density approximation (LDA). The flatness condition implies that the 
trapping potential has to vary slowly compared to the typical local de Broglie 
wavelength of the particles. As a result, LDA can be seen as a WKB or a 
semiclassical approximation. 

If we apply the condition for the validity of the WKB approximation |116) 
along a certain direction € {x, y, x} in the trap, we obtain 



where the trapping potential is given by V cx (x 1 y, z) = (muj 2 x 2 + mKj/ + 
mu; 2 z 2 )/2, the harmonic oscillator length in the direction of r» yields U = 
(hjmijji) 1 ! 2 , and where the Fermi energy eF CT (r) = /z. 2 (67r 2 n CT (r)) 2 / 3 /2m sets 
locally the typical kinetic energy scale for the atoms. For the early MIT ex- 
periments, we have ui x = u> y — 2ir x 115 s _1 in the radial direction, and 
ui z = 2ir x 22.8 s _1 in the axial direction. Moreover, the total number of 
trapped atoms is about 10 7 , leading to central densities on the order of 10 18 
m -3 [147] . Plugging in the numbers, we see that the MIT experiments indeed 
easily satisfy the condition for the local-density approximation in both the ra- 
dial and the axial direction. Only near the edge of the cloud, where the atomic 
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densities become very small, the approximation breaks down. However, for the 
early Rice experiments, we have that ui x = ui y = 2ir x 350 s _1 and uj z = 2ir x 7.2 
s _1 , while their total number of particles is about a 10 5 [25]. As a result, a 
question mark arises over the validity of the local-density approximation, when 



they study very large imbalances. Namely, then the left-hand side of Eq. (35 1 
is not much smaller than unity anymore for the minority density in the steep 
radial direction. 

However, in the rest of this section we focus more on the MIT experiments, 
namely we assume that it is valid to use the LDA. The approximation then con- 
veniently allows us to absorb the trapping potential in the chemical potentials 
and use locally the homogeneous theory that was developed in the previous sec- 
tion. Moreover, we assume that we have conveniently rcscalcd the coordinate r, 
by a factor Q/uJi with Q = (uixUjyijJz) 1 ^ 3 , so that the external potential becomes 
spherically symmetric. As a result, we have for the spherically symmetric trap 
in rescaled units that (J, a (r) = — V cx (r) — \i a — mij 2 r 2 /2, which leads for 
the local average chemical potential to fi(r) = fi— y cx (r) = \l — mui 2 r 2 /2, while 
(half) the difference remains constant, i.e. h(r) = h. Moreover, we can calculate 
the total particle number N a with spin a in the trap by 



/drn CT [(A)(r);T,/v(r)], (36) 



where the local particle densities n CT [(A)(r); T, Mo-(r)] are given by Eq. (28). 
Since at position r in the trap the homogeneous phase is realized that corre- 
sponds to T/fi(r) and h/fi(r), we can already predict what the trapped config- 
urations look like by considering only the homogeneous phase diagram of Fig. 
T2^a). Namely, we may draw an arrow in the homogeneous diagram that pre- 
cisely follows those ratios T//i(r) and h/fj,(r) that are encountered in the trap 
|146j . We put the tail of the arrow at T//i(0) and h//i(0) corresponding to the 
center of the trap. For increasing r or decreasing //(r), both ratios increase and 
the arrow consequently moves to the right and to the above in the diagram, so 
that we end up in the normal phase. Only when both h = T = 0, we stay in 
the origin of the diagram and superfluidity is encountered throughout the trap. 
For h = and T ^ 0, the arrow moves upward along the vertical axis and we 
encounter a second-order superfluid to normal transition in the trap. 

However, we are interested in the polarized case, so that we consider h > 0. 
Then, the arrows can follow six different types of paths, corresponding to six 
trapped configurations, as shown in Fig. |13| The first case is when the whole trap 
is in the normal phase as follows from arrow (a). The corresponding schematic 



representation of the trapped configuration is depicted on the right in Fig. 13 
In this representation, we schematically show which phases are encountered as 
a function of the radius r in the trap. The thick circle then represents the edge 
of the cloud. The second case is indicated by arrow (b), where the cloud has a 
BCS superfluid core. For increasing r, we now first encounter a crossover to the 
gapless Sarma regime, after which there is a second-order phase transition to 
the normal phase. The third case, arrow (c), leads to a gapless superfluid core, 
which is separated from the normal outer region by a second-order transition. 
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The fourth case, arrow (d) , has a fully gapped superfluid core, surrounded by a 
Sarma shell, which is separated from the normal edge by a first-order transition. 
The fifth case, arrow (e), has a gapless superfluid core, and the transition to the 
normal state in the trap is of first order. The sixth case, arrow (f), is similar to 
the fifth case, only the core is now fully gapped. 

In the next section, we present the universal phase diagram for the strongly 
interacting Fermi mixture in a trap. Although we have specified 6 trapped 
'phases', some of them only differed by the presence of a crossover. Since a 
crossover connects different physical regimes in a smooth manner, these regimes 
are strictly speaking not different thermodynamic phases due to the absence of 
a true phase transition. In that sense the Sarma phase and the BCS phase are 
two different limits of the same polarized superfluid phase at nonzero tempera- 
tures. In the BCS limit, the polarization is caused by thermal excitation of the 
gapped quasiparticle branches, whereas in the Sarma limit the gapless branch 
in principle does not need temperature to be filled. If we choose to classify the 
trapped phases solely in terms of the true phase transition that occurs in the 
trap, only three different trapped phases remain. First of all there is the normal 
phase, given by arrow (a), for which there is no transition in the trap. The sec- 
ond phase has a second-order transition in the trap, such that it encompasses 
both arrows (b) and (c). Note that the phases (b) and (c) have in common 
that there is always a region of Sarma superfluidity present. Therefore, we call 
the phase with a second-order transition in the trap also the Sarma phase. The 
third trapped phase has a first-order transition in the trap and thus encompasses 
arrows (d), (e) and (f). Due to the presence of a first-order interface, we also 
call it the phase-separated phase. The three different trapped phases are also 
illustrated in Fig. [8j 

3.6. Phase diagram in a trap 

The main results of the mean-field calculations in the trap are presented 
in Fig. [14] Here, we show the universal phase diagram of a trapped Fermi 
gas in the unitarity limit as a function of temperature and polarization |85] . 
This phase diagram is universal in the sense that it does not depend on the 



total number of fermions or the trap geometry. Fig. 14 reveals that there is 
a tricritical point for the trapped Fermi mixture, which was also the case for 
the homogeneous gas [84 ) I146|. 1162] . As explained in the previous section, the 
normal phase means that the gas is in its normal state throughout the trap. In 
the Sarma phase the Fermi gas has a shell structure, in which the core of the 
trapped gas is superfluid, whereas the outer region is normal. Furthermore, the 
normal-to-superfluid transition as a function of the position in the trap is of 
second order. Since the superfluid order parameter (A) vanishes continuously 
at the transition, we have for nonzero polarizations always a region in the trap 
where |(A)|(r) is so small that it results in a gapless superfluid with negative 
quasiparticle excitation energies, as first studied by Sarma [51]. Since |(A)|(r) 
increases towards the center of the trap, it is also possible that the superfluid 
becomes gapped in the center of the trap. This leads to a gapped BCS superfluid 
core with a gapless Sarma superfluid and normal shell surrounding it [148] . 
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Figure 14: Universal phase diagram of a trapped imbalanced Fermi gas in the unitarity 
limit. The polarization P is given by (N+ — N—)/(N+ + N-), where N± designates the 
number of fermions in each hyperfine state of the Fermi gas. The temperature T is scaled 
with the critical temperature T c of the balanced trapped Fermi gas and the polarization is 
scaled with the critical polarization at zero temperature P c . Using mean-field theory, we have 
that /cbT c = 0.66/i with fi the average chemical potential at the center of the trap, while 
P c = 0.998. 



Finally, in the phase-separated region of the phase diagram, the superfluid core 
and the normal shell of the gas are separated by a first-order transition as a 
function of position, which implies that (A)(r) goes discontinuously to zero at 
a certain equipotential surface in the trap. 

Fig. [M] allows for a natural explanation of the qualitative differences in the 
early observations by the experimental groups at MIT [147) and Rice university 
[2"5] . More precisely, we could argue that the initial experiments of Zwierlein et 
at have observed the smooth transition from the normal phase to the Sarma 
phase at a large polarization |147j . implying that these early experiments have 
been performed above the temperature of the tricritical point. Moreover, we 
then suggest that the experiments of Partridge et al. have been performed in 
the temperature regime below the tricritical point, since these experiments ap- 
pear to see the transition between a non-phase-separated and a phase-separated 
superfluid phase at small polarization [25 . Note that this explanation based on 
mean-field theory can only qualitatively account for the major initial discrep- 
ancies between the Rice and the MIT experiments, but fails on a quantitative 
level to describe the strongly interacting experiments accurately. Morefore, as 
mentioned in Section |3.I[ to fully understand the early results by the Rice ex- 
periment non-equilibrium effects must be taken into account [73 1151] . 

Now, we turn to the question how we actually obtained the phase diagram 
of Fig. [14] We first determined the line between the normal and the two super- 



42 



fluid trapped phases. This is achieved by chosing two central chemical potentials 
jLto-(O) and determining in the center of the trap the temperature at which the 
BCS order parameter vanishes. To this end we use the homogeneous theory of 
Section [3^41 where inspection of the BCS thermodynamic potential revealed that 
the vanishing of the order parameter can occur continuously or discontinuously, 
i.e., by a second-order or a first-order phase transition. If the transition is of 
second order in the center of the trap, we go from the trapped normal phase 
to the trapped Sarma phase. If the transition is of first order in the center of 
the trap, we go from the trapped normal phase to the trapped phase-separated 
phase. At the tricritical point in the center of the trap these two different kinds 
of transitions merge. This procedure thus gives us two central chemical poten- 
tials and a critical temperature of the gas, either above or below the tricritical 
point. We now still have to calculate the corresponding total atom number and 
polarization with the use of Eq. ( 36 ) in order to find a point in the phase di- 
agram of Fig. [l4j that lies on the line between the normal phase and the two 
superfluid trapped phases. Note that it is convenient that the phase diagram is 
universal, so that the diagram is independent of the total atom number. This 
feature saves us an additional iteration procedure compared to the case when 
it is necessary to calculate the phase diagram for a certain specific total atom 
number. 

So far, we looked primarily at the center of the trap, but the tricritical 
condition can also be satisfied at a certain equipotential surface outside the 
center of the trap. This gives us a point on the Sarma-to-phase-separation 
line [85]. To see this, consider a point on this line and raise the temperature 
slightly. This changes the tricritical transition outside the center of the trap 
into a second-order transition slightly closer to the center of the trap, which 
means that the gas is in the Sarma phase. In a similar way, a slightly lower 
temperature leads to a first-order transition as a function of position in the 
trap, i.e. to the phase-separated phase. So the procedure now gives us by 
construction a tricritical temperature and two chemical potentials at a certain 
radius r c3 outside of the trap center. To calculate the corresponding total atom 
number and polarization we use again Eq. (36), where we note that for r < r c ^ 
the gas is superfluid, while for r > r C 3 the densities are normal. 



3. 7. Mass-imbalanced case 

So far, we have explored the mean-field phase diagram of the strongly- 
interacting Fermi g function of temperature and spin polarziation. Most 
recently, experiments have indicated that the physical consequences of yet an- 
other parameter can be explored, namely that of a mass imbalance between the 
two fermionic components. A very promising mixture in this respect consists of 
6 Li and 40 K, which has a mass ratio of 6.7. Several accessible Feshbach reso- 
nances are identified in the mixture jSU [95] , while both species have also been 
simultaneously cooled into the degenerate regime [551127]. In the unitarity limit, 
the size of the Cooper pairs is comparable to the average interparticle distance 
and the pairing is a many-body effect. The mass imbalance has a profound 
effect on the pairing, because it affects the way in which the two Fermi spheres 
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overlap. Theoretical studies at zero temperature using mean-field theory [163J 
and Monte-Carlo simulations |164| considered the BEC-BCS crossover in this 
system. We show next that for the large mass ratio of the 6 Li- 40 K mixture, 
the phase diagram in the unitarity limit is even richer at nonzero temperatures 
than for the mass-balanced case. Similar to the solely spin-imbalanced case 
is the presence of phase separation |163l 11651 155] , which can occur due to the 
mismatch of the Fermi surfaces. Also similar is that gapless Sarma superflu- 
idity is unstable at zero temperature S5J, while there is a predicted crossover 
to the Sarma phase at nonzero temperatures [89] . However, the most striking 
difference is the presence of a Lifshitz point in the phase diagram [9"9"] . 

At a Lifshitz point the transition to the superfluid phase undergoes a change 
of character. Rather than preferring a homogeneous order parameter, the sys- 
tem now becomes an inhomogeneous superfluid. This exotic possibility was 
early investigated for the weakly interacting mass-balanced case by Larkin and 
Ovchinnikov (LO), who considered a superfluid with a single standing- wave or- 
der parameter 103 . This is energetically more favorable than the plane-wave 
case studied by Fulde and Ferrell (FF) [104J. Since the LO phase results in 
periodic modulations of the particle densities, it is a special kind of supersolid. 
[94] I108j . The FF and LO phases have intrigued the physics community for 
many decades, but so far remained elusive in experiments with atomic Fermi 
mixtures. Typically, Lifshitz points are predicted at weak interactions where 
the critical temperatures are very low [1001 11011 11021 [75] . However, the phase 
diagram of the 6 Li- 40 K mixture contains both a Lifshitz and a tricritical point 



in the unitarity limit, as shown in Fig. 15 This is in sharp contrast to the 
mass-balanced case, where mean-field theory only leads to a tricritical point at 
unitarity, which is in agreement with experimental observations [72] . 

As we discussed in Section [3T4] , the critical properties of the superfluid tran- 
sition in a fermionic mixture are determined by the effective Landau theory for 
the superfluid order parameter A(r). Although we consider no external poten- 
tial, the order parameter may still vary in space due to a spontaneous breaking 
of translational symmetry. Close to the continuous superfluid transition, we 
expand the effective grand-canonical thermodynamic potential as 

0[A] = |dr| 7L |VA| 2 + aL |A| 2 + ^|A| 4 |, (37) 

where the challenge is to express the expansion parameters in terms of the tem- 
perature T and the atomic chemical potentials /i± with the upper (lower) sign 
referring to the light 6 Li (heavy 40 K) atoms. In this section, we use again mean- 
field theory to achieve this goal. A phase transition has occurred when the global 
minimum of f2 is located at a nonzero order parameter (A(r)), which describes 
a condensate of pairs. When 7l is positive, the pairs have a positive effec- 
tive mass and their center-of-mass state of lowest energy is at zero momentum. 
Then, we can consider a homogeneous pairing field A, for which a second-order 
transition occurs at a critical temperature T c determined by ccl(Tc) = 0, for 



which the expression was given in the mass-balanced case by Eq. (33). For the 



mass-imbalanced case, the same expression holds, but now the (average) kinetic 
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Figure 15: Universal phase diagram for the homogeneous 6 Li- 40 K mixture in the unitarity 
limit as a function of temperature T and polarization p. The temperature is scaled with the 
reduced Fermi temperature k^Tp = ep = h 2 (3ir 2 n) 2 ^ /2m, where m is twice the reduced 
mass and n is the total particle density. The phase diagram is calculated with mean-field 
theory. For a majority of light 6 Li atoms there is a tricritical point (TCP), at which the 
normal state (N), the gapless supcrfiuid Sarma phase (S), and the forbidden region (FR) meet 
each other. For a majority of heavy 40 K atoms there is a Lifshitz point (LP), where there 
is an instability towards supersolidity (SS). The size of the supersolid stability region is not 
calculated within our theory and the dashed lines are therefore only guides to the eye. The 
dash-dotted line indicates the crossover from gapless Sarma to gapped BCS superfluidity. 



energy is given by et = fi. 2 k 2 /2m with m = 2m+m_/ (rn + + m_) and the Fermi 
distributions become / CT (k) = l/[ e ^( e ^k-AO _j_ with e CT .k = fi 2 k 2 /2m CT . 

A continuous transition only occurs when the vanishing minimum at A = 
is the global minimum, which is not necessarily the case. As discussed in Section 



3.2 the expansion of may contain higher powers of |A| that have negative 



coefficients, leading to a first-order transition with a jump in the order parameter 
when f2[0] = f2[(A)]. Second-order behavior turns into first-order behavior when 
/3l becomes negative, so that the temperature T C 3 at the tricritical point (TCP) 
is determined by a\,(T c z) = and /3l(?c3) = 0. For the mass-imbalanced 
case, the expression for /?l is given by Eq. (34) with the same expressions for 
m and /o-(k) as were just introduced for a^. Another possibility is that not 
but rather 7l goes to zero. This leads to a Lifshitz point (LP), which is 
thus determined by «l(2l) = and 7l(7l) = 0. Since the effective mass of 
the Cooper pairs becomes negative below the Lifshitz point, it is energetically 
favorable for them to have kinetic energy and form a superfluid at nonzero 
momentum. This can be established in many ways, namely through a standing 
wave [103 or a more complicated superposition of plane waves [12. 93, 102, 80 . 
Due to the variety of possibilities it is hard to predict which lattice structure is 
most favorable. However, the fact that they all emerge from the Lifhitz point 
facilitates the experimental search for supersolidity in the 6 Li- 40 K mixture. The 
supersolid phase could be observed using Bragg spectroscopy. 

The mean-field BCS thermodynamic potential density for the 6 Li- 40 K mix- 
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ture is given by Eq. (27 1 with £k the average kinetic energy, m twice the re- 
duced mass, and the Bogoliubov quasiparticle dispersions given by Hui^^ — 
fiwk — cr[2h — e + .k + e_,k]/2 with fiwk = \J (ek — m) 2 + |A| 2 . We can apply 
the mentioned critical conditions to wbcs an d obtain the mean-field phase dia- 
gram. Although the BCS potential neglects fluctuations in the order parameter, 
it is expected that these fluctuation effects only result in quantitative correc- 
tions. We have seen this already for the strongly interacting experiments in 
the mass-balanced case, where the mean-field diagram seems to be qualitatively 
correct |72] . The coefficients determining the second-order phase transition and 
the tricritical point are readily calculated from the thermodynamic potential 
as a L = <9w B cs[0; av]/<9|A| 2 and /3 L = <9 2 cj B cs[0; /v]/<9 2 |A| 2 , resulting in the 



mentioned expressions. The obtained phase diagram [99] is shown in Fig. 15 
where the polarization is defined as p — (n + — + n_) with n + (ft.-) 

corresponding to the Li (K) density. These particle densities are determined 
by n a = — <?Wbcs[(A); n a \/dn„. Therefore, the polarization is discontinuous 
simultaneously with the order parameter, which gives rise to a forbidden region 
(FR) below the tricritical point. 



From Fig. 15 we see that the mean-field phase diagram also contains a 
Lifshitz point. It can be calculated from the noninteracting Green's function 
(i.e. the propagator) for the Cooper pairs G&{iu} n ,\s), which is discussed in 



more detail in Section 4.2 From Eq. (37), we see that 7 L can be interpreted 



as an inverse effective mass of the Cooper pairs, while — ai, corresponds to a 
Cooper-pair chemical potential, and /3l can be seen as a Cooper-pair interaction. 
In the normal state, the Cooper-pair propagat or is given by hG^} (iui n , k) 



1/T — hH{iuj ni k), which is derived in Section 4.2 leading to Eq. (55 1. Here, 
T 2b = Anah 2 /m with a the diverging scattering length, while S is the expression 
for the so-called ladder diagram, shown in Fig. [24] It is given by 

w ^ f dk ' J 1 l-/+(k')-/_(k-k') | 

ha(iu n , k = / — — - <^ — — + — — \ . (38 

J (2tt) j [2e(k') ihuj n - e +M , - e_, k -k' + 2^ J 

We have that the mean-field expression for «l is equal to — /z,G^ 1 (0, 0), while 
7l = — dhG^{Q 1 0)/9k 2 = at the Lifshitz point. By applying this condition 
we find that the Lifshitz point occurs for a majority of heavy 40 K atoms, which 
can be qualitatively explained in the following way. The ideal case for pairing is 
an equal amount of particles with the same mass, so that the two Fermi surfaces 
fully overlap. In the case of a mass difference, the Fermi energy of the lighter 
particles is larger than that of the heavier particles, while the Fermi wavevectors 
stay the same. If consequently the number of heavy particles is increased, then 
the difference in Fermi energies is reduced, enhancing the tendency for pairing. 
Since now the Fermi wavevectors are different, the pairing is expected to occur 
at nonzero momentum. If the number of light particles is increased, then the 
difference in Fermi energies is further enhanced resulting ultimately in phase 
separation. There is a critical mass ratio for which the phase diagram changes 
its form from having two tricritical points (like for the mass-balanced case) to 



a tricritical and a Lifshitz point (like in Fig. 15). At unitarity, it is given by 



4G 



r = 4.2 according to mean-field theory. At this multicritical point we thus have 
that «l = = 7l = 0. A continuously tunable mass difference might be 
achievable with two different atomic species in an optical lattice. 

What precisely happens below the Lifshitz point is an open question for fur- 



ther research [12]. In Fig. 15 we have sketched a simple scenario, where there 



is a second-order transition from the normal to the supersolid phase, for which 
the condition is G^ 1 (0, k s ) = with k s the wavevector of the supersolid. How- 
ever, this scenario might not be sufficient since the transition could in principle 
also be of first order, where the realized supersolid periodicity can contain more 
than one wavevector [93] . The transition from supersolidity to the homogeneous 
superfluid phase is also expected to be of first order. The calculation for the 
stability regions of all possible supersolid lattices is beyond the scope of this 
review, so that the dashed lines in Fig. 15 are merely guides to the eye. A more 
elaborate discussion on various possible supersolid phases is given in Ref. [12] . 

To include the presence of an external trapping potential, the local-density 
approach can be used for the (local) superfluid phases, like the BCS phase or 
the Sarma phase. This gives rise to shell structures for the phases in the trap, 
just like for the mass-balanced case in Section 
case such an approach was pursued in Ref. [ 



3.6 For the mass-imbalanced 



66] . However, to treat (non- 



local) supersolid phases in the trap we need theories that go beyond LDA. 
The Bogoliubov-de Gennes approach, which will be treated in Section |6.2[ is a 
well-established approach to treat non-local superfluid phases and it has been 
applied to the mass-imbalanced case at zero temperature in Ref. [167] . Another 
approach that not only goes beyond LDA, but also beyond mean-field theory is 
the Monte-Carlo approach, which was applied to few-particle systems by Blumc 
[168J. This method is restricted so far to zero temperature. 



4. Diagrammatic approaches 

In Section [3] we theoretically studied the strongly interacting Fermi gas with 
a population imbalance using the mean-field BCS thermodynamic potential. In 
this section, we start from the microscopic Hamiltonian that describes the in- 
teracting atomic Fermi gas and derive diagrammatic approaches that go beyond 
mean-field theory. We mainly use the functional formalism, in which the cen- 
tral object is the fermionic action that belongs to the microscopic Hamiltonian. 
The diagrammatic approaches in the functional formalism give rise to a versatile 
theoretical toolbox. This allows us to understand the main shortcomings of the 
mean-field approximation, which already was adequate to qualitatively describe 
several aspects of the experiments in the strongly interacting regime. In this 
section, we try to improve the theoretical calculations so that also quantitative 
agreement with experiments is reached. 

The starting point is the Hamiltonian for interacting fermions in two dif- 
ferent spin states, which we label with a spin index a = ±. As mentioned 
before, for the atomic population-imbalanced Fermi gas the spin label repre- 
sents two different internal hyperfine states of the atomic ground state. Since 
the resulting hyperfine space is only two dimensional, it is often also referred 
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to as a (pseudo)spin-l/2 space. We consider utracold atoms that interact via 
s-wave interactions, which is the dominant scattering mechanism at the very 
low energies and momenta of interest, as explained in Section [2] There, it was 
also shown that s-wave scattering only occurs between fermions in different spin 
states. Considering the point interaction V(r — r') = VoS(r — r') to model the 
short range interactions between the atoms, we can write down the following 
second-quantized Hamiltonian in the grand-canonical ensemble, 



H = J rfr <^( r ){- 

dr y o ^(r)^(r)0_(r)0 + (r) . (39) 



Here, the creation operator ^J-(r) (annihilation operator </>cr(r)) creates (anni- 
hilates) a particle with mass m a in spin state a at position r, where \i a is 
the corresponding chemical potential. In this section, we discuss interacting 
particles with identical mass so m+ = m_ = m. However, the two chemical 
potentials can be different in order to account for a difference in the population 
of the two spin states. 



To the second-quantized Hamiltonian of Eq. ( 39 1 corresponds a microscopic 
action £ [</>], which is obtained from a derivation that generalizes Feynman's path 
integral approach to many-body quantum physics, see e.g. Ref. |169j . The big 
advantage of the functional or path-integral approach is that the operators <P\{t) 
and <%.(r) make place for fields 0*(r, r) and <^ CT (r, r), which are often easier to 
deal with in calculations. The variable r is then an (imaginary) time variable, 
that can be used to incorporate either the dynamics, or, as in our case, the equi- 
librium quantum statistics of the system. The position vector is given by r. The 
fields are ordinary complex numbers in the case of bosons and anticommuting 
Grassmann numbers in the case of fermions. More precisely, the fields 4> a { T i r ) 
correspond to eigenvalues of the annihilation operators, whose eigenstates are 
called coherent states. Working with the eigenvalues of the operators is often 
convenient, because it is typically easier to manipulate numbers than operators. 
Following the derivation in Ref. |30j , we have that the microscopic action in the 
functional formalism becomes 

S M = -E / dTdT ' / drdr ' €(^r)^i(r,r;r',r')^(r',r') 

rH/3 r 

+VoJ o drj dr^(r,r)^*(r,r)^_(r,r)^ + (r,r), (40) 

where we have introduced the so-called non-interacting (inverse) Green's func- 
tion HGq] t1 given by 

hG^(r, r; r', r') = _ /ft JL - ^ - *( r - r ')S(r - r'). (41) 
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The Fourier transform of the inverse Green's function yields 

hGo^(iu) n ,k) = ihw n - e k + (42) 

where = h 2 \<L 2 /2m is the kinetic energy, and w„ are the odd fermionic Mat- 
subara frequencies, as seen for example in Refs. [1571 150] . 

4--1- Hubbard- Stratonovich transformation 

As mentioned in Section [3j BCS theory can be physically interpreted as the 
Bose-Einstein condensation of Cooper pairs. This follows from the observation 
that the order parameter for the BCS transition is proportional to the expec- 
tation value of the pair annihilation operator (0_(r)0+(r)), which means that 
the Cooper pairs are in a coherent state analogous to a condensate of point- 
like bosons. For the transition to the paired condensate, we require that the 
two-body interaction potential is attractive, because otherwise the formation of 
pairs would not be energetically favorable. From now on, we therefore consider 
an attractive potential. Note that this does not necessarily mean that the cor- 
responding scattering length a of the interaction is negative. In Section [2] we 
namely saw that attractive potentials can also give rise to a positive scattering 
length when there is a two-body bound state present in the potential. To intro- 
duce the BCS order parameter exactly into the many-body theory, we use the 
so-called Hubbard-Stratonovich (HS) transformation |170| I17R Il72j . 

To perform this transformation, we start from the expression of the partition 
sum Z in the functional formalism, whose derivation can be found for example 
in Ref. We have that 

Z = [V(f> e - s M/ h , (43) 



which represents a functional integral over all fermion fields <j> and <jf that are 
antiperiodic on the imaginary time interval [0, In the functional integral of 



Eq. (43), we can insert the following identity 



1 = y^expj^ drydr[A*(r,r) -Vb^(r, r)^*(r,r)] 

x 7" 1 [A(r,r) - 7o0_(r,r)^ + (r,r)] /fcj , (44) 

which represents a functional integral over all bosonic fields A and A* that 
are periodic on the imaginary time interval [0, fi/3]. Here, Af is a normalization 
constant that cancels the outcome of the Gaussian functional integral. Since 
this constant does not depend on any thermodynamic variable, it turns out to 
be irrelevant for our purposes. For notational conveniency, we will absorb this 



numerical constant in the measure of Eq. (44). Moreover, from Eq. (44) it is 



also clear that on average the complex pairing field A(t, r) satisfies 

(A(r,r)) =V (4>-(r,r) ( f> + (T,r)). (45) 
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More details and a more general discussion about the HS transformation is given 
inRefs. [17211371]. 



By combining Eqs. (43 1 and (44), we obtain 

Z = [v<j)VA e - s ^ A ]M 



(46) 



where the new action is given by 



S[0,A] 



Fill 



dr / dr 



A(r,r) 
V 



(47) 



rH/3 p 

h J dr dr' J dr dr' &( T , r) • G^ s (t, r; r', r') ■ *(t', r'). 



Here, we have introduced ^(t, r) = [<jft (r, r) , (f)_ (r, r)] and the 2x2 inverse 
Green's function matrix GgJ s in spin space, which is given by 

G B i s (r,r;r',r') = G " 1 (r, r; r', r') - £ B cs(r, r; r', r'), (48) 

with the noninteracting part 



Go 1 (r,r;r',rO = 
and the self-energy part 

^S B cs(T,r;r',r') = 



Go^(r,r;r',r') 




A(r,r) 
A*(r,r) 





-Go^(T,r;T',r') 
<5(r-r')(5(r-r'). 



(49) 



(50) 



The minus sign in the matrix element G^\ 2 01 Eq. (49) comes from reversing 
the order of the two corresponding fermionic Grassmann fields in Eq. ( 47 1 . The 
pairing field A couples to two fermionic creation fields in Eq. (47), showing 
that a pair can decay into two fermions of opposite spin, while A* couples to 
two fermionic annihilation fields showing that two fermions can form a Cooper 
pair. The action of Eq. (47) can be physically interpreted as an interacting 
Bose-Fermi mixture. We see that by performing the exact HS transformation, 
we have eliminated the fourth-order interaction term in Eq. ( 40 1 , so that the 
resulting action of Eq. (47) is only quadratic in the fermion fields. However, 
this goes at the cost of introducing an additional functional integral over the 
pairing field A. 

Since the action of Eq. ( 47 ) is quadratic in the fermionic fields, the corre- 
sponding functional integral is Gaussian and can be performed exactly. Using 
the standard formula for sermonic Gaussian functional integrals, see for example 
Rcf. 130 , we obtain 



Z 



with the action given by 



S[A] 



dr / dr 



A(r,r) 
V 



ffRpog(-G^ s )]. 



(51) 



(52) 
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Here, the trace is to be taken over spin space, imaginary time and position space, 
or, after Fourier transformation, over spin space, frequencies and momenta. 
We thus see that with the HS transformation, we have been able to derive 



two actions that are equivalent to Eq. ( 40 ) , and that thus exactly introduce the 



order parameter field into the theory. Since none of the three actions are exactly 
solvable, we have to invoke approximations in order to get actual results. How- 
ever, having three actions at our disposal greatly enhances the scope of possible 
approximation schemes. Moreover, it often depends on the physical quantity of 
interest which action is most suitable as a starting point. For example, if we 
are interested in the fermionic self-energy in the normal state, we might start 



from Eq. (40 1, while if we are interested in the properties of the Cooper pairs, 



we could start from Eq. ( 52 ) 



By expanding the logarithm in Eq. (52 1 in powers of A, we obtain an infinite 



series. As just mentioned, this prohibits an exact solution to the problem, and 



approximations have to be made in order to proceed. The part in Eq. (52 1 that 



is independent of A gives rise to the partition sum of the ideal Fermi gas in 



Eq. (51). By performing the mean- field, or saddle-point, approximation, the 
full path integral over the bosonic field A(r, r) is simply approximated by the 
value of the integrand associated with the global minimum of the action. This 



approximation results in the BCS grand potential of Eq. (27). For a detailed 
derivation, see for example Ref. 30 . 

4-2. Pair fluctuations 

Although mean-field BCS theory provides an elegant description of supercon- 
ducting electrons in metals or superfluidity of weakly interacting atomic Fermi 
gases, it has some major shortcomings in the strongly interacting regime. In 



the normal state, when (A) = 0, the BCS grand-potential density of Eq. (27) 
reduces to the grand-potential density of the ideal gas. This is because the 
BCS theory only takes the interaction of (quasi)particles with the condensate of 
Cooper pairs non-pertubatively into account, but not the interaction between 
the (quasi) particles themselves. Once the condensate has vanished in the normal 
state, we thus have that wbcs(0; T, n a ) = Wj g (T, with 

c ig [A;l> CT ] = * + (53) 

k,<7 

As a result, BCS mean-field theory does not lead to an accurate description of 
the normal state in the unitarity limit, because the particles are then actually 
very strongly interacting, as follows from the diverging scattering length. 

Another issue is that BCS theory describes only formation of pairs with zero 
center-of-mass momentum, i.e., it describes only a condensate of Cooper pairs. 
The critical temperature following from BCS theory then corresponds to the 
typical temperature at which the pairs are broken. This is in particular a bad 
approximation in the BEC limit of the crossover that was discussed in Section 



2.5 In the BEC limit the superfluidity is namely lost at the critical temperature 
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due to thermal occupation of non-zero momentum states by the (tighly-bound) 
pairs, that are still intact. 

Both problems of the BCS theory can already be overcome to a significant ex- 
tent by considering also the Gaussian fluctuations of the order parameter around 
its mean-field value. Since Gaussian functional integrals can be performed ana- 
lytically, the contribution of these pair fluctuations can thus be studied exactly, 
resulting in a theory that was pioneered by Nozieres and Schmitt-Rink [129 . 
Diagrammatically this theory corresponds to summing the so-called ladder dia- 
gram, or ring diagram contributions to the grand potential. Summarizing their 
analysis, we start in the normal state and expand the logarithm in Eq. (52 1 to 
second order in the pairing field, as also more elaborately explained in Ref. [30 . 
The Cooper pair Green's function G A X is then defined through this quadratic 
part of the Cooper pair action, and is found to be given by 

HG^^r-r'y) = -b(r - r')S(r - r') 

+ ^G , + (r, r; r', r')G 0) _(r, r; r', r'), (54) 
or, after Fourier transformation, by 

l-/+(e k ')-/-(ek-k') 1 \ (55) 



AittPa V \ —ihuj n + £k' + £k-k' — 2/i 2et' 



where we eliminated the microscopic interaction strength Vg in favour of the 
two-body scattering length a using Eq. ( 23 ) . 

The resulting functional integral after substituting the quadratic part of the 
Cooper pair action in Eq. (51) can be evaluated exactly by using the formula 
for bosonic Gaussian functional integrals, sec e.g. Ref. [169 . The partition sum 
Z up to quadratic order is found to be given by a part that is independent 
of the non-condensed Cooper pairs, and a part coming from the Gaussian 
integration over the Cooper pair propagator, Z A , which can thus be physically 
interpreted to describe a noninteracting gas of Cooper pairs. The latter is given 
by 

Z A = e - T ^- G ±^W\ (56) 

while Zq is given by the ideal gas partition sum for fermions 129, 30J . By using 
Z = Z Zj\ and the thermodynamic relation Z = e~P Vu «°( T < IJ "'') with cj gc (T, /v) 
the grand-canonical thermodynamic potential density, we thus conclude that 

w gc (T,/i CT ) = w ig (T,/v) + u A (T,n„), (57) 

where the contribution due to the Cooper pairs uj&(T, is given by 

wa(T,/v) = ^^log[-G A 1 ( l ^„,k)] (58) 

n.k 

= ^E/ dw &Mlm(log[-G A V+,k)]), 
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with 6(e) = l/(e' 3£ — 1) the Bose distribution function and uj + = ui 
rj is infinitesimally small and positive 



where 



In the second step of Eq. (58), we have 
used contour integration in order to write the sum over Matsubara frequencies 
as an integral along the real axis |129j . This is convenient because the result- 
ing frequency integral converges very fast. Moreover, the imaginary part of 
G^ 1 (o; + ,k) along the real axis can be determined analytically and yields 



Im(G A V + ,k))= 



8(hio + 2zt - e k /2)m 2 



CT=±1 



47T/?|k|/l 4 



(59) 
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where 6{x) is the Heaviside step function. The real part of G^ l (uj + , k) can then 
be determined by a Kramers-Kronig transform. Note that to this end, we have 
to subtract the analytically known high-frequency tail of the imaginary part, 
which is given by (haj + 2/i — e^/2) 1 ^ 2 m 3 ^ 2 /inh 3 . Only then, the remaining part 
decays fast enough to be transformed. 

For the spin-balanced Fermi gas, the Nozieres-Schmitt-Rink approach has 
been applied succesfully to the whole crossover from the BCS regime to the BEC 
regime [130] . The theory gives rise to a critical temperature curve that smoothly 
interpolates between the result from BCS theory in the weakly-interacting limit 
and the result for Bose-Einstein condensation of noninteracting bosons in the 
extremely strongly interacting limit, where a two-body bound state is present. 
The critical condition for superfluidity within this theory is that G^ 1 (0, 0) = 0, 
which can physically be interpreted as the chemical potential for the Cooper 
pairs being equal to zero. This thus corresponds to the well-known criterium 
for condensation of bosons. By comparing Eq. (55) with Eq. (33), we see that 
the critical condition has not changed from mean-field theory, resulting in the 
unitarity limit in k^T c = 0.66/i c for the critical temperature. This is rather large 
compared to most recent measurements by Ku et al., that give rise to ks,T c = 
0.40/x c [T75] . or to Monte-Carlo results, that predict k B T c = 0.31/i c [174] . But 
since the equation of state has changed due to the inclusion of the non-condensed 
Cooper pairs in Eq. (57), we find by calculating the total particle density at the 
critical temperature with n(T, /z) = — duj gc (T, /z)/ctyz that the NSR theory leads 
to T c = 0.23Tp, which is a significant decrease from mean-field theory. The 
latter prediction is rather close to the results from various more sophisticated 
strong-coupling theories, that typically predict critical temperatures in the range 
T c = 0.15±0.05T F [T7J[I7J[TIB1[I771[T7H], while the most recent experimental 
value is T c = 0.167 ± 0.013T F (173] . 

Using Eqs. (57) and (58) we can calculate the pressure of a strongly inter- 
In the unitarity limit, 



-U! 



acting Fermi gas, given by p s (T, /z CT ) 
this pressure was measured for a balanced Fermi gas by Nascimbene et al. |74j , 
Horikoshi et al. |179j , and by Ku et al. [173] . We show the results of Nascimbene 
et al. and of the prediction by the NSR theory in Fig. 16 where the agreement 
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Figure 16: (Color online) Equation of state for the normal phase of a strongly interacting 
balanced Fermi gas at unitarity in the grand-canonical ensemble. The pressure p g of the gas is 
calculated as a function of the inverse fugacity e~ ^ M with the renormalization group approach 
of Section |5.5| (full line) and the Nozieres-Schmitt-Rink approach (dashed line). The pressure 
of the ideal two-component Fermi gas is given by pi g . The triangles are the experimental 
results of Nascimbene et al. 1741 . 



is seen to be very good. We note that the experimental results of Nascimbene 
et al. and Ku et al. agree reasonably well with each other for the temperature 



range of Fig. 16 while the results of Horikoshi et al. |179j do not match the 
high-temperature virial expansion |180j . At lower temperatures the results of 
Nascimbene et al. and Ku et al. differ quite significantly, where the results 
of Ku et al. seem to agree better with most recent Monte-Carlo calculations 
|181j . It is possible to generalize the NSR theory to the superfluid state by 
making in Eq. (47) the substitution A(r, r) = Ao + A'(t, r), expanding Eq. (52) 
up to second order in the fluctuations A', and performing the resulting Gaus- 
sian functional integral [122, 182, 183J . In the superfluid state, the NSR theory 
gives rise to good agreement with earlier experiments 184, 185J, while the most 
recent experimental data of Ku et al. agrees best with diagrammatic Montc- 
Carlo calculations [181] . A self-consistent fluctuation theory that takes fcrmionic 
self-energy effects into account also leads to rather good agreement with these 
recent experiments 1861 1175] , From this discussion, we conclude that current 
high-precision measurements allow for a very sensitive test of many-body the- 
ories. Although the NSR theory is not exact, and therefore is expected to give 
rise to deviations from experiments in the unitarity regime, it is seen from Fig. 
|16| to give a big step forward compared to mean-field theory. 

The fact that the NSR theory of pair fluctuations is not perfect is most 
clearly seen at nonzero polarizations, where it gives rise to unphysical results. 
Namely, for small imbalances the NSR theory predicts near the critical temper- 
ature a negative polarization p = (n + — n_)/(n + + rt_) for a positive chemical 
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potential difference — /it_), which corresponds to a compressibility matrix 
— d 2 uj gc (T, Ha) /d n a d n a i that is not positive definite |146j . It is a quite unsat- 
isfactory situation that the NSR theory, which gives such good results for the 
balanced Fermi gas, already gives unphysical results for very small population 
imbalances. In Section [5.5| we discuss a way to overcome this problem. 



4-. 3. Fermionic self- energy 

So far, we have seen that the pair fluctuation effects in a strongly-interacting 
Fermi gas lead to an accurate description of the thermodynamics for a balanced 
Fermi gas, but fail to describe to imbalanced Fermi gas. An extreme example of 
an imbalanced Fermi gas is a single spin-down particle interacting with a Fermi 
sea of spin-up particles, also called the Fermi polaron [15ULI187] . The self-energy 
of the Fermi polaron has been calculated with various methods including 
Monte-Carlo calculations [TSD1 HHH1 Q2S] , variational methods [T571 [HM HUT] , 
the many-body T matrix in the ladder approximation |192| . and RG methods 
|193j . The various theories are in good agreement with each other, and with 
experiment [73 11194] , leading to about KE_ — —0.6fi+ at zero temperature. We 
note that this large self-energy effect is not incorporated in the pair-fluctuation 
theory of the previous section. 

The goal of the present section is to discuss in more detail the strongly- 
interacting normal state of the imbalanced Fermi gas at zero temperature as 
described by the equation of state from the pioneering Monte-Carlo study of 



Lobo et al. |150j . The results of their calculations are shown in Fig. 17 The 
circles present the Monte-Carlo data, while the full line is the best fit to this 
data. The dashed line represents the following Ansatz by Lobo et al. for the 
total energy density e of the interacting system, 

3 3 m 3 

e = -n + €F+ + - — n_eF- — -Aep+n- 
5 5 m* 5 

= 5 rl +e F+ (l + ^ 5/3 -^), (60) 

which was constructed to treat the strongly polarized regime. This Ansatz 
physically describes a few spin-down atoms in a Fermi sea of spin-up atoms, 
where the minority atoms have acquired an effective mass of m* = 1.04m and 
a self-energy of 7i£_ = — 3Aep+/5 with A = 0.97. Both effects are the result 
of the strong attractive interactions with the spin-up Fermi sea, which remains 



unaffected. The Ansatz of Eq. (60) is less suitable to describe the balanced 
case, x — 1, where it causes for example the two chemical potentials \i a to be 
unequal. 

To have an Ansatz that is suitable for all polarizations, we may write the 
equation of state from the beginning in a form that is symmetric upon inter- 
change of n + and ra_. This reflects the fact that it does not matter for the 
energy density which of the two spin states is the majority state. We try a 
self-energy of the form 

hX„(n + ,n_) = — „ )1/3a , (61) 
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Figure 17: Equation of state for the zero-temperature normal state of the spin-imbalanced 
Fermi mixture with unitary interactions. The Monte Carlo results by Lobo et al. |150] are 
shown by circles, while x = n—/n+. The full line represents the best fit of Lobo et al. to their 
data. The dashed line gives the result for considering only the self-energy of the minority 
particles. Upon taking also the self-energy of the majority particles into account through a 
symmetric ansatz, the dashed-dotted line is obtained. The latter is hard to discern from the 
best fit. The energy density e is scaled with the ideal gas energy for the spin- up particles 
<Ei g+ = 3e P+ n + /5 with e FtT = h 2 (6ir 2 n a ) 2 ^ s /2m. 



where we still have to specify the power a. If we use the self-energy of Eq. (61 1 



to obtain a total energy density similar to Eq. ( 60 1 , we find 



3 3 

-n + e F+ + -n_e F _ 
5 5 



3A(6tt 2 ) 2 / 3 ^ 2 



r 5/3 



.4 



10™ 

X 



(1 + a; Q ) 1 /3" 



(62) 



which agrees for small x with Eq. ( 60 ) . If we wish to obtain agreement with 
the full Monte-Carlo equation of state, we may use A = 1.01 and a = 2, which 
results in the dashed-dotted line shown in Fig. [17] It is hardly discernible from 
the full line that gives the best fit of Lobo et al. The value of A is calculated in 
the next section, while there does not seem to be a direct microscopic argument 
why a should be 2. It can thus be interpreted as a single fit parameter. 



4-. 3.1. The extremely imbalanced case 

We have just seen that the equation of state for the normal state at zero 
temperature is to a large extent determined by self-energy effects of the minority 
(spin-down) atoms in the Fermi sea of majority (spin-up) atoms. In this section, 
we calculate the self-energy of a single spin-down atom in a Fermi sea of spin-up 
atoms analytically in the so-called ladder or many-body T matrix approximation 
as was done by Combescot et al. [192] . The effect of this self-energy fi.E (T (/i cr ) is 
to renormalize the chemical potential according to ^(^a) = /v — hY, a (Ha)- We 
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may let the renormalized chemical potentials fj! a consequently enter the grand 
potential density of the normal gas at zero temperature w n in the following way 

uj n (fi a ) = w ig (/^.(T = 0,/v)) (63) 




The single-particle limit for spin state a is located right at the transition from a 
nonzero particle density to a zero-particle density. The latter takes place when 
the renormalized chemical potential fj,' a goes to zero, i.e. when fi a = fiE^/v). 
To be able to solve this equation, we need to calculate the self-energy K£, a (fjL a ) 
in the extremely imbalanced limit. 

To perform the self-energy calculation for a single spin-down atom in a Fermi 
sea of spin-up particles. The many-body T matrix is treated in the so-called 
ladder approximation, which is for example described in Ref. |30) . For unitary 
interactions and at zero temperature, the many-body T matrix in the extremely 
imbalanced case is seen to have a particularly simple expression, because the 
Fermi distributions become step functions and there is no spin-down particle 
density, resulting in 

where 9(x) is the Heaviside step function. The T matrix gives rise to a self- 
energy for the minority particle, which at zero momentum and frequency is 
given by [H2] 

KE-fo) = / dwT(w,k')Go, + (w J k') 

k' ~ io ° 

= ^E T ( e k'-M+' k ')^+-e k <), ( 65 ) 

k' 

where the substitution hio — > — fi + in the T matrix comes from a contour 
integration over the frequency uj. Then, to fulfill the initial assumption of zero 
down particles, we still need to solve n'_ = 0, which leads to /i_ = — 0.607/i+ = 
— 3A/i + /5 with A = 1.01. This result means that for a chemical potential 
lower than \i_ = — 0.6/i + there are no spin-down particles, whereas for a higher 
chemical potential there is a nonzero density of spin-down particles. Of course, 
the results would be the same for a single spin-up particle in a spin-down sea. 
There have been several Monte-Carlo calculations for a single fermion with spin 
a in a Fermi sea of particles with spin —a. The results are A = 0.97 |150j . 
A = 0.99 [US] and A = 1.03 [HE], which all agree remarkably well with the 
simple ladder calculation. In the following, we will keep using A = 1.01. 

The fact that the T matrix calculation in the ladder approximation agrees so 
well with Monte-Carlo results can be regarded to some extent as a coincidence 
(190] . Namely, we could also try to perform a more selfconsistent calculation, 
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Figure 18: Diagrammatic representation of the so-called bubble sum, leading to a screened 
two-body transition matrix T sc . 



in which the fermionic self-energy effects are also taken into account in the 
expression for the many-body T matrix. Such a selfconsistent T matrix approach 
was pioneered by Haussmann for the balanced Fermi gas |186j . In order to take 
fermionic self-energy effects into account, we could replace /i_ by fi'_ = in 
Eq. (64). However, then we find that /i_ = ft£- = —0.9/^+, which does 
not agree with the Monte-Carlo results. To proceed further we could try to 
incorporate another effect, namely that of particle-hole excitations that tend 
to screen the interaction, as is well known for electrons in a metal |154j . We 
may account for this screening by replacing the two-body transition matrix T 2b 
with an effective transition matrix that includes the particle-hole excitations 
through the so-called bubble sum, as is depicted in Fig. [18] We then have 
that l/T sc (ui, k) = 1/T 2b — hH(ui, k), where the expression for the polarization 
bubble diagram is given by, see e.g. [30] . 



MI(w,k) 



/ + (k + k')-/_(k') 



k' 



2h' - e k+k > 



(66) 



If we now replace 1/T 2b with 1/T SC (0, 0) and /z_ by fi'_ = in Eq. (64), then we 
find \X- = US- = —0.5{i+, which is again quite close to the mentioned Monte- 
Carlo results of /i_ = — 0.6^ + . Here, we have neglected the momentum and 
frequency dependence of the screened interaction. It actually turns out that 
nonzero external momenta reduce the effect of screening, where we note that 
a slight reduction of the screening effect would bring the present calculation 



even closer to the Monte-Carlo results and recent measurements (194] . We may 
conclude from this discussion that in order to accurately describe the normal 
state of the imbalanced Fermi gas, in particular the fermionic self-energy and 
the effect of screening seem to be of importance. 

We then try to generalize the fermionic self-energy effect on the polaron 
to the full equation of state by using an Ansatz for the normal state that has 
the correct extremely imbalanced limits in it, but actually leads to excellent 
results with Monte-Carlo calculations for all polarizations. This is achieved by 
introducing the following renormalized chemical potentials 



3 A- 
5 



M-CT 



(67) 



We can now determine readily the renormalized chemical potentials n' a as a 
function of the microscopic ones /i CT , and consequently use n a = —duj n / djjL a to 
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Figure 19: (Color online) Equations of state for the zero-temperature normal state of the 
B Li- 4t) K mixture (full line) and the mass-balanced mixture (dashed line) in the unitarity limit. 
The Monte Carlo results for the 6 Li- 40 K mixture by Gezerlis et al. |195) are shown by circles, 
while the Monte Carlo results for the mass-balanced mixture by Lobo et al. I150| are shown 
by squares. The energy density e is scaled with the ideal gas result e; g = 3(ep+n+ + ep_n_)/5 
with t-ptr = ^ 2 (67r 2 n a -) 2/ ' 3 /2m [T . The deviations from one thus show the strong interaction 
effects. 



determine the densities. The two extremely imbalanced solutions to Eq. (67 1 
are such that the chemical potential of the majority species (— a) is not renor- 
malized, = /Lt_ CT , while the renormalized chemical potential of the minority 
species a is zero, n' a = fi a + 3A[i- a /5 = 0. It is easy to invert the quadratic re- 
lations in Eq. (67) to obtain f/ a (fj,a-), which motivates the choice for a quadratic 
Ansatz. However, a better validation of the approach is obtained by comparing 
with the Monte-Carlo equation of state. We determine the energy density using 



e = uj n + n+n + + //_n_, for which the result is shown in Fig. 17 Here, both the 



result of the present procedure (dashed line) and the Monte-Carlo data of Lobo 
et al. |150j (squares) are shown, giving very good agreement. Moreover, we 
also applied the same procedure to the case of the strongly-interacting 6 Li- 40 K 
mixture [99], where we compared with the Monte-Carlo data of Gezerlis et al. 
[195J. Again, the agreement is excellent. 

The approach is also expected to be appropriate at small nonzero temper- 
atures, since the coefficient A is not expected to be strongly temperature de- 
pendent at low temperatures in the normal phase. This is because the Fermi 
sea of majority particles remains unaffected for T < Tp, where the Fermi tem- 
perature typically represents a large energy scale in the system. It is given by 
fcjjTjr = £f = ft 2 (37r 2 7z) 2//3 /2m with n = n + + n_ the total atomic density. We 
can check this assumption by comparing with the Monte-Carlo calculation at 
nonzero temperature of Burovski et al. In this calculation, the chemical po- 
tential at the critical temperature in the unitarity limit was determined to be 
fi(T c = 0.15Tp) = 0.49£f |174j . We can compare the last result with the direct 
extension of the present approach to nonzero temperature, which follows from 
using ui n (X 1 , fi a ) = o;bcs(0;T,/4). By applying n a = 5w n (0.15T F , fi a )/dn a , we 



59 



0.5 1 0.5 

a) r/R+ b) r /R + 



Figure 20: Observation of phase separation in a trapped Fermi mixture by Shin et al. |72| . 
The measurement was performed at their lowest temperature of about T = 0.02Tp + (0) and 
a total polarization of P = 0.44. Their data are shown by the noisy curves, a) Calculated 
density profile for the spin-up particles n+ (full line), the spin-down particles n_ (dashed 
line), and their difference (dashcd-dotted line). The densities are scaled by the central density 
of the spin-up particles and R+ = (2fi+ /mui 2 ) 1 ^ 2 . b) Exactly the same, but now for the 
column densities, which means that an additional integration along the y axis is performed. 
Even though the calculations are performed at zero temperature, the agreement is very good. 



obtain /Lt(0.15Tp) = 0.53eF in quite good agreement with the Monte-Carlo cal- 
culation. However, we note that both results differ somewhat from the most 
recent measurements for the chemical potential as a function of temperature, 
which had a maximum at /Lt(0.17Tp) = 0.42ep [173] . 

4-3.2. Density profiles 

Having obtained the grand potential density for a homogeneous Fermi gas in 
the strongly interacting normal state, we can now use it to also study the normal 
trapped gas by means of the local-density approximation, which was discussed 
in Section [3~5| Namely, since we have that the local chemical potential for spin 
state a is given by n<j{r) = fi a — V cx (r), where [i a is the chemical potential 
in the trap center, we immediately obtain the particle density at each point in 
the trap by using n a (r) — — du) n (T, /v(r))/9jU CT . In the case that the gas cloud 
does not have a superfluid core, this procedure gives rise to the complete density 
profile for the trapped normal phase in the unitarity limit. 

However, if the spin imbalance is not too large, then the core of the trapped 
unitary Fermi mixture is superfluid. In Section [3j we discussed the mean- 
field BCS theory for a superfluid Fermi gas. At zero temperature, BCS the- 
ory predicts for the equal-density superfluid equation of state at unitarity that 
fjb = (1+/?bcs) £ f with /?bcs = —0.41 |136j . If Gaussian fluctuations are included 
in the superfluid state, we have that /3nsr = —0.60 196 , while Monte-Carlo 
calculations give /3m c = —0.58 [197] . Measurements give results ranging from 
/3 cxp = -0.49 to -0.68 PHEZlUHElUS]. In the followin g we use the Monte-Carlo 
value of /3mc = —0.58. From the equation of state, we can immediately calculate 
the energy density of the superfluid. Namely, from the relation /i = de s f/dn, we 
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obtain upon integration that e s f = 3(1 + /3mc)£f^/5. As a result, we find for 
the grand potential density in the superfluid phase at zero temperature that 

2 

w s f(/v) = e si - f^n = --(1 + /3 MC )e F n (68) 

5 

2(2m)3/V M 

— -(1 + PMCj 



15tt 2 /i 3 V 1 + # 



When we have that oj s f(/i cr ) = u; n (/z CT ,0), then a first-order phase transition be- 
tween the equal density superfluid and the polarized normal phase occurs. From 



Eqs. (63 1, (67) and (68 1, it follows that the critical difference in the chemical 
potentials is given by h c = 0.93/1. At this critical value, the (local) polarization 
jumps from zero in the superfluid phase to p c = (n_|_ — n._)/(n+ + n_) = 0.38 
in the normal phase as found by Lobo et al. [150] . Note that this result is very 
different from the mean-field result of Section [3~4j where we obtained p c = 0.93 
for the critical polarization. The difference shows the crucial quantitative effect 
of the strong interactions in the normal state. In the trap, the equation of state 
of Lobo et al. in combination with the local-density approximation gives rise to 
P c = 0.78 150 , as opposed to P c = 0.998 that follows from mean-field theory, 
see Fig. [14] 

By comparing lu s { with w n everywhere in the trap, we can determine which 
of the two phases is locally most favorable and calculate the particle densities 
from the corresponding equation of state |199j . For the superfluid phase, the 
local particle densities are obtained from /z(r) = (1 + ^Mc)ep( r )- We are now 
thus in the position to determine the density profiles for the trapped phase- 
separated state, which has a superfluid core and a first-order transition to the 
normal state as a function of position in the trap. This discontinuous transition 
has been observed in the density profiles measured by Shin et al. |72j , as shown 



in Fig. 20 In panel (a), the density profiles are shown for the spin- up particles, 
the spin-down particles and their difference. The noisy curves are the data, 
which compare well with the calculated lines. The experimental profiles were 
obtained at a temperature of about T = 0.02Tp+(O), where Tp+(0) is the local 
Fermi temperature of the spin-up particles in the center of the trap. The total 
polarization was about P = 0.44 with a total number of spin-up particles given 
by N + = 6 • 10 6 . Moreover, we specified the trapping potential V ex (p, z) = 



m(oj\(x 2 + y 2 ) + lo 2 z 2 )/2 for the MIT experiments in Section 3.5 From Fig. 
20^ a), we see that the measured densities are equal in the superfluid core, but 
upon going outwards from the center of the trap there is a sudden rise in the 
density difference, indicating a first-order transition. In the figure, we used the 
radial coordinate r given by w 2 r 2 = ljj\(x 2 + y 2 ) + uj 2 z z 2 with u = (w^^) 1 / 3 . 
The calculated profiles are shown by the full line for n + (r), by the dashed line 
for n_(r), and by the dashed-dotted line for their difference. The densities 
are all scaled by the central density of the spin- up particles n+(0). The radial 
distances are scaled in the figure by the radial cloud size of the spin-up atoms, 
given by R + = (2fi + /mui 2 ) 1 / 2 . 
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In Fig. 20 [b), also the column densities shown, which follow from 

poo 

n ca (x,z) = / dy n a (x,y,z). (69) 



The column densities are less noisy, because they are directly probed by the 
in-situ imaging that is performed in the experiment, while the densities have 
to be reconstructed. The radial coordinate for the column density is now given 
by w 2 r 2 = lo\x 2 + uj 2 z 2 , while the plotted column densities are all scaled by 
the central column density of the spin-up particles n c+ (0, 0). In general, the 
agreement between the experimental curves and the theoretical curves is very 
good. The small differences could be due to the fact that the calculation was 
performed at zero temperature, while the experiment is at very small nonzero 
temperature. Moreover, the local-density approximation is expected to fail close 
to the interface, because a true jump in the order parameter would cost an 
infinite amount of gradient energy. If gradient terms are taken into account, 
then a smooth order-parameter profile is calculated, as we will see in more 
detail in Section [6] 

Note that throughout the calculation in this section, we have assumed that 
more exotic superfluid phases, like the FF and LO phases mentioned in the 
introduction, or induced p-wave superfluidity 92 , do not play a role at zero 
temperature in the unitarity limit. To put it differently, we have assumed that 



the mean-field phase diagrams of Figs. 12 and 14 are topologically correct. Seen 



the good agreement between experiments and theory, this assumption seems to 
be justified so far. 



5. The renormalization-group approach 

In the following section, we discuss in some detail the renormalization-group 
(RG) approach to interacting quantum gases in order to dicuss in a different 
way the importance of fermionic self-energy and screening effects. Renormaliza- 
tion group theory is not only a powerful technique for studying quantitatively 
strongly-interacting gases, but it also gives an elegant conceptual framework 
for understanding many-body physics and phase transitions in general. The 
latter comes about because we are often interested in determining the physics 
of a many-body system at the macroscopic level, i.e., at long wavelengths or 
at low momenta. As a result, we then need to determine the effect of micro- 
scopic degrees of freedom with high momenta to arrive at an effective quantum 
field theory for the long-wavelength physics. This is achieved elegantly by the 
momentum-space renormalization group approach as formulated originally by 
Wilson 200, 201]. In this section, we study his approach because of its simplic- 
ity both conceptually and technically. However, nowadays many formulations 
of the renormalization-group equations are available |202l 12031 12041 1205j , having 
each their advantages and disadvantages. The renormalization group procedure 
and related ideas have been also applied to the unitary Fermi gas, where early 
examples include the use of the e expansion [206] . the 1/N expansion [1771 1207] . 
and the functional renormalization group approach 208, 209J. 
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The goal of Wilsonian renormalization is to construct a transformation that 
maps an action, characterized by a certain set of coupling constants, to a new 
action for longer wavelengths, where the values of the coupling constants have 
changed. This is achieved by performing two steps. First, an integration over 
high- momentum degrees of freedom is carried out, and the effect of this inte- 
gration is absorbed in the coupling constants of the action that are therefore 
said to flow. Second, a re-scaling of all momenta and fields is performed to 
bring the relevant momenta for the action back to their original domain. By 
repeating these two steps over and over again, it is possible to arrive at highly 
nonperturbative approximations to the exact effective action. 

At a continuous phase transition, the correlation length diverges which im- 
plies that critical fluctuations are dominant at each length scale and that the 
system becomes scale invariant. This critical behavior is elegantly captured by 
the renormalization-group approach, in which a critical system is described by 
a fixed point of the above two-step transformation. By studying the proper- 
ties of these fixed points, it is possible to obtain accurate predictions for the 
critical exponents that characterize the nonanalytic behavior of various thermo- 
dynamic quantities near the critical point. Moreover, the renormalization-group 
approach also explains universality, which is the observation that very different 
microscopic actions give rise to exactly the same critical exponents. It turns out 
that these different microscopic actions actually flow to the same fixed point, 
which is to a large extent solely determined by the dimensionality and the 
symmetries of the underlying theory. As a result, critical phenomena can be 
categorized in classes of models that share the same critical behavior. 

In this review, however, we do not wish to calculate universal critical expo- 
nents, since these have already been determined with great accuracy for the XY 
universality class, which is the class of the transition to the superfluid state in 
three spatial dimensions. Rather, we want to calculate quantities like the critical 
temperature of a unitary interacting Fermi gases, that have only been measured 
in recent years. Such quantities are not universal from a renormalization group 
point of view, and depend on the microscopic details of the Hamiltonian. There- 
fore, in this section we use the renormalization group only as a nonperturbative 
method to iteratively solve a many-body problem, rather than as a map be- 
tween actions with the same high-momentum cutoff from which critical scaling 
relations can be derived. As a result, the renormalization step in which the 
fields and coupling constants are rescaled is not particularly useful for our goals 
and we leave it away. 

5.1. The first step of Wilsonian renormalization 

Thus, we only perform the first step of the renormalization method, which 
is to evaluate the Feynman diagrams that renormalize the coupling constants 
of interest, while keeping the integration over the internal momenta restricted 
to the considered high-momentum shell. Only one-loop and tree diagrams con- 
tribute to the flow, because we consider the width of the momentum shell to be 
infinitcsimally small, and each loop introduces an additional factor proportional 
to the infinitesimal width. Although the one-loop structure of the infinitesimal 
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Wilsonian RG is exact, it does not mean that it is easy to also obtain exact 
results, since this would require the consideration of an infinite number of cou- 
pling constants. Since the latter is usually not possible in practice, the RG 
distinguishes between the relevance of the coupling constants, so that a small 
set of them may already lead to useful non-perturbative results. 

To start the derivation of the Wilsonian renormalization group equations, we 
consider the action in momentum space, S = So + Si , where the noninteracting 
part is given by 

So[4>] =-Y\ #,n,k^oi(* w n> k )6r,n ) k, ( 70 ) 



with the noninteracting Green's function of Eq. (41 1. The interacting part is 
given by 

V 

Sl[<f>] = Y ^+,n',k' < / , -,m-n',q-k'^-,m-n,q-k^+,r 1 ,k , (71) 

k,k',q 

ro, n ,m 

where the notation with the prime, namely Vq, indicates that the interaction 
strength is a renormalizing quantity in our treatment. In the same way, we 
will use the notation \J a for the renormalizing chemical potentials. To arrive 
at the one-loop corrections to this action arising from the integration over fluc- 
tuations in an infinitesimal high-momentum shell, it is particularly convenient 
to use the following procedure, where all one-loop corrections can be obtained 
by performing a single Gaussian functional integral. Namely, by splitting the 



functional integral of Eq. ( 43 1 into a part that contains the integration over 



low-momentum modes </>„ k with small wavenumbers |k| < A, and a part that 
contains the integration over high-momentum modes <fi* k with large wavenum- 
bers in an infinitesimal shell A — dA < |k| < A, we obtain 

< e -s [<l> < ]/h f T>lh > p -s [4> > ]/h p ~s 1 [<t,<,<t>>]/h 
2?0< e -S[* < ]/» 

J P0>ex P | Y, *Sk' ' G-\iu n/ ,k') ■ *>A 
j P^e-^l^expj-Trllog^-G- 1 )]}, (72) 



where we introduced the notation &^J k , = 



b +*n<M>>< 



for the fast 



modes in the high-momentum shell. In the second line of Eq. ( 72 1 , we moved 
all terms that only depend on long- wavelength modes (</>^ k ) to the first in- 
tegral, while in the exponent of the second integral on the third line we only 
kept the terms that are quadratic in the short- wavelength modes (</>^, k , ). We 
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note that other terms, like of linear order or higher orders in </>>, would lead 
after expansion of the exponent and integration over (/>> either to tree diagrams 
or to higher-loop Feynman diagrams. The tree diagrams that are allowed by 
momentum conservation give rise to effective three-body interactions or higher, 
which we ignore in the present discussion. Moreover, we argued that higher-loop 
diagrams vanish because we work in the lim it o f infinitesimal shell width. The 
inverse Green's function matrix G 



where G is the Fourier transform of Eq. (491, while we have for S> that 



_1 of Eq. (72) is given by GZ 1 



h 2 /3V 



n.k 



< * 1, 

+ ,n,k^— , — n, — k 



^-,-n,-k^+,n.k 
^+,n,k^+,n,k 



(73) 



which follows directly from Eq. (71), Eq. (72) and assuming that the interaction 
stays local during the RG flow. 

In the last step of Eq. ( 72 ) , we have performed an exact integration over 



the fast modes, resulting in the well-known expression for a Gaussian functional 
integral. The trace is to be taken over the 2x2 spin space structure of G _1 , 
over all Matsubara frequencies, and over the momenta in the considered high- 
momentum shell. To expand Tr[log(— G^ 1 )] in powers of (/><, we use 



log(-G: 



log(-G - 1 )+log(l~G S > ). 



(74) 



For example, expanding the logarithm to first order in S>, we find, after taking 
the trace over the 2x2 Nambu space, the contributions proportional to (/><*<^<. 
Namely, we have that 



k , n 



k.n 



n,k^0 2-J 
k' ,n' 

n ^%M^A^ 

k.n 



Go,+ (k>') 



(75) 



where in the second line we first summed the Green's function over the fcrmionic 
Matsubara frequencies resulting in the Fermi distribution / ff (k) = l/{exp[/3(ek— 
li'v)] + 1} (see e.g. Ref |157] for more details), after which we converted the sum 
over k' into an integral. Since we integrate over a shell of infinitesimal width, 
the outcome is simply the value of the integrand at the momentum hA times 



the width of the shell dA. Note that Eq. (75 ) is valid for integration from A to 
A + dA, i.e., integrating out fluctuations from lower to higher momenta, where 
the cutoff A sets the scale up to which fluctuations have already been integated 
out. As a result, we would have an additional minus sign for integration of 
fluctations from higher to lower momenta, since then dA becomes negative. 



Eq. ( 75 ) can be interpreted as a self-energy contribution, which can be ab- 



sorbed in the 'renormalized' chemical potential of the minority particles, //_• 
We thus find that 

d f / a = -^V>f- a (A)dA, (76) 
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Figure 21: Feynman diagrams renormalizing a) the chemical potentials and b) the interatomic 
interaction. 



which is a differential equation for \j! a . The change in the chemical potential 
is seen to be proportional to the (renormalized) interaction strength V ' and 
the (renormalized) average density of particles, as expected for a self-energy 
contribution. By expanding Tr[log(— G^ 1 )] to second order in S>, we can also 
find the renormalization of the interaction strength, resulting ultimately in 



-Vn 



A*_ 

2tt 2 



1-/+(A)-/_(A) /+(A)-/_(A) 
2(e A - /i') 2h' 



dA, 



(77) 



with fjf = + //_)/2 and h = (// + - n'_)/2. The nrst term in E Q- < 77 l 
corresponds to the so-called ladder diagram and describes the scattering be- 
tween particles. The second term corresponds to the so-called bubble diagram 
and describes screening of the interaction by particle-hole excitations. Also 
note that due to the coupling of the differential equations for p! a and Vq 1 , we 
automatically generate an infinite number of Feynman diagrams, showing the 
nonperturbative nature of the RG. The expressions in Eqs. (76) and (77) are 
diagrammatically represented in Fig. |21| The only difference here compared to 
the usual Feynman rules is that the internal momenta in the one-loop diagrams 
are restricted to stay in the infinitesimal high-momentum shell. 



5.2. Renormalization theory for fermions 

From the previous discussion, we conclude that if we consider the interaction 
vertex to remain frequency and momentum independent during the RG proce- 
dure, and if we only take into account the renormalization of the three coupling 
constants fi' a and Vq, we find the following coupled RG equations [193J 

dA ~ 2tt 2 
d^ = A 2 
dA 2tt< 

with the Fermi distribution f a (A) = l/{exp[/3(eA — ^f a )] + 1}. To be able to 
solve these coupled differential equations numerically, we need to specify the 
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MO 



2h' 



(78) 
(79) 



66 



initial conditions. We start the RG flow at a high- momentum cutoff HAq, where 
the chemical potential is not yet renormalized, so that we have fjf a (Ao) = [i a 
with \x a the (not renormalized or bare) chemical potential, which is a given 
thermodynamic variable when working in the grand-canonical ensemble. As 
an initial condition for the interaction, we use Vq'(Aq) = Air 2 ati 2 / m(ir — 2aA ) 
with a the s-wave scattering length. We encountered this expression already 
before in Eq. (24). There, this expression for the two-body case was found 
to be the exact relation between the microscopic (not renormalized or bare) 
interaction parameter Vq and the effective interaction strength, namely the two- 
body transition matrix T 2b with scattering length a. When we take the two- 
body limit of Eq. (78 1, i.e. f a — > and [i! a — > 0, we find the differential form 
of Eq. (23). This implies that the two-body limit of the RG equations gives 
after integration rise to a renormalized interaction strength of T 2h — Airah 2 jm. 
The initial condition for the interaction is thus chosen to make sure that we 
automatically incorporate the correct two-body scattering length into the RG 
theory. Moreover, this initial condition also ensures that at the end of the RG 
calculation no dependence on the arbitrary high-momentum cut-off fiAo remains. 
The results depend on the thermodynamic variables fj,„, T and the scattering 
length a. In the unitarity limit, when a — > oo, we have that Vq (Ao) = 
—mA /27r 2 h 2 , so that then also the dependence on a disappears. 

Although we have specified the initial conditions, the outcome of the differ- 
ential equations is not yet completely fixed. Because the differential equations 
are coupled, this final outcome depends on the specific way we perform the in- 
termediate steps, i.e. it depends on the way in which we flow. Ultimately, we 
want to take the effect of the fluctuations in all momentum shells into account, 
and to this end we have to specify which momentum shells to consider first and 
which to consider last. If we would be able to calculate all (infinitely many) 
coupling constants generated by the RG procedure exactly, then the order in 
which we integrate out the fluctuations would not matter anymore, because all 
effects of the fluctuations in one momentum shell on the other shells would be 
fully included. However, for our current simplest nontrivial RG equations, this 
is not the case. If the RG is formulated for classical systems or bosonic systems, 
then usually the most relevant excitations of lowest energy have a momentum 
near zero. Therefore, the RG equations are then taken to start at the high- 
momentum cutoff Ao, and the effects of an increasing amount of fluctuations 
on the low-momentum modes is considered by using for example A(Z) = A e~ l . 
If we insert A(/) and dA(l)/dl into Eqs. (78) and (79), we get a set of coupled 
differential equations in I which we can numerically evolve from I = Q to I = oo. 
Note that then an additional minus sign is needed, because we flow from high 
to low momenta. 

For a fermionic system the relevant excitations of lowest energy have a mo- 
mentum near the Fermi energy, which is therefore a natural endpoint for a RG 
flow |2I0j . In an imbalanced system, actually three Fermi levels are relevant, 
namely one for each species and the average Fermi level. Moreover, these Fermi 
levels, given by the renormalized chemical potentials n' a {l), are shifting during 
the flow. If we would want to flow automatically to the renormalized average 
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Fermi level, we could use the following procedure. First, we start at a high 
momentum cutoff HAq and flow to an intermediate cutoff hA' to integrate out 
high-momentum (two-body) physics. Then, we start integrating out the rest of 
the momentum shells approximately symmetrically with respect to the flowing 
average Fermi level to treat the many-body physics. This is achieved by using 

urn 



and by 



a+(o = K - 
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2m//' (I) 

h 2 



2m/i 



2m//(Z) 



(80) 



(81) 



h 2 V h 2 

Note that, as desired, we have that A + (Z) starts at A' and automatically flows 



from above to y/2mfi'(oo)/h, whereas A_(Z) starts at and automatically flows 
from below to WZmf7(po)/h. This behavior is illustrated in Fig. 22 By substi- 
tuting A + (l), A- (I) and their derivatives in Eqs. (78 1 and (79), we obtain a set 



of coupled differential equations in Z that again can be solved numerically. Note 
that this procedure leads to different results than from using A(Z) = A n e~ l . A 
similar procedure can also be used for the imbalanced case to flow automatically 
to either one of the Fermi levels /^.(Z), or even to both Fermi levels simultane- 
ously. 

As a somewhat technical side comment we note that the two functions hA a (l) 
are not exactly symmetric with respect to y/2mfi'(l). This can be a problem 
when there are true poles in the RG equations, as happens for example to 



Eq. (78) at zero temperature. The RG equation for the interaction can then be 



interpreted as the calculation of the principal value for a shifting pole. In order 
to get finite results the divergent parts on each side of the pole have to exactly 
cancel each other, which can be somewhat cumbersome to achieve, since also 
the derivative of the shifting pole plays a role. This can be seen from Eqs. (78), 



(80) and (81) and noting that dA±/dl enter the differential equation on both 
sides of the pole. However, for our present purpose this problem does not play a 
major role, since we are interested primarily in nonzero temperatures and then 



Eq. (78) has a well-defined right-hand side for all A(Z). 



Unfortunately, close to the critical temperature, Eq. (78 1 gives rise to another 
problem. Namely, when the Fermi mixture is critical, the inverse many-body 
vertex Vq flows to zero according to the Thouless criterion |1291 I211| and 



the chemical potentials in Eq. (79) consequently diverge, which is unphysical. 
The most simple way to solve this problem and calculate critical properties 
more realistically, is to take the momentum dependence of the renormalizcd 
interaction vertex into account. Namely, considering the ladder diagram in 



Fig. (21 ), we note that in general this diagram depends on the external center- 
of-mass momentum q and w m of Eq. (71). As a result, this diagram leads in 
general to the contribution 



dZ(q 2 
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dk l-/ + (e k )-/_(e q _ k ) 
dA ( 27r ) 3 -itWm + £k + e q _ k - 2ju' 
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Figure 22: a) Position of the momentum shells (dashed lines) and flow of the chemical 
potentials (solid lines) for a) the strongly interacting extremely imbalanced case, b) the weakly 
interacting balanced case, and c) the strongly interacting imbalanced case. 



where during integration over high-momentum modes both k and q — k have 
to remain in the infinitesimal shell of width dA. We note that the expression on 
the right-hand side of Eq. (82) depends only on q 2 , rather than on q. 



We thus have that the renormalization procedure can lead to a momentum- 
dependent interaction, even though the microscopic interaction is to a very good 
approximation momentum independent. If we consider only ladder diagrams, 
which is a commonly used approximation for dilute quantum gases because 
ladder diagrams do not vanish in the two-body limit, then the renormalizcd 
interaction will depend in first instance only on q and uj m . In order to simplify 
things even further, we only take the center-of-mass momentum into account by 
expanding the (inverse) interaction in the following way: V' 1 = V ' — Z' a q 2 



The flow equation for the additional coupling constant Z' q 

by [na 



dZ' q 



ddZ(q 2 ,0) 



"0 

is then obtained 



(83) 
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5.3. The extremely imbalanced case 

To apply the discussed RG formalism, we consider first one spin-down parti- 
cle in a Fermi sea of spin-up particles at zero temperature in the unitarity limit, 
also called the Fermi polaron, which was discussed in Section 4.3.1 When we 



calculate the self-energy of the polaron with the present renormalization group 
approach, the RG equations are simplified, because /-(A) can be set to zero 
and thus is not renormalized. In Section 4.3.1 it was mentioned that the 



contour integration to determine the one-loop Matsubara sum in the diagram 
of Fig 
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a), leads to the frequency substitution ifmj r , 



//, in Eq. (82) 



Using this and accounting for the external momentum dependence with the 
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coupling Z' , we obtain |193| 
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(85) 
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Eq. (84) is similar to Eq. (78), where the main difference is the denominator of 



the ladder diagram contribution caused by the just described substitution that 
takes into account the frequency-dependence of the interaction. We note that 
in the above equations we have not considered the effect of external momenta 
and frequencies on the bubble diagram contribution, which is an approximation. 



Eq. ( 85 ) is similar to Eq. ( 79 ) , but the denominator has changed because we 



include the effect of the momentum-dependent interaction Vq in the self-energy 
contribution. As a result, we have that the interaction in Eq. (75 1 would now be 



given by Vk + ic' rather than Vq. This also means that the resulting renormaliza- 
tion contribution to the self-energy depends on k, giving rise to a renormaliza- 
tion of the effective mass of the fermions. However, in our present treatment we 
only consider the renormalization of the chemical potential, i.e., the momentum- 
independent part of the self-energy, and we put k = 0. Eq. (86) follows from 



Eqs. (82) and (83), where we set /_ to zero and perform the mentioned fre- 
quency substitution. The initial conditions for the above coupled differential 
equations are Vq~ 1 (A ) = —mA /2ir 2 h 2 , //_(0) = \i- and Z' q ~ (0) = 0, since 
the interaction starts out as being momentum independent. 

Note that for the case to our interest we have that //_ is initially negative 
and increases during the flow due to the strong attractive interactions. The 
quantum phase transition from a zero density to a nonzero density of spin-down 
particles occurs for the initial value /i_ that at the end of the flow precisely 



leads to /i'_(oo) = 0. We find fj,. 



-0.52^ + , where we used A(Z) = A e to 



flow to the final value of the renormalizcd minority Fermi surface. Note that we 
thus again have to include the additional minus sign for going from high to low 
momenta. This initial value of fjf_ is therefore also the self-energy of a strongly 
interacting spin-down particle in a sea of spin- up particles |192j . Our value 
for the self-energy is somewhat less than the results obtained by Monte-Carlo 
calculations [150| I188[ I189j . giving rise to /i_ = — 0.6//+, while the experiment 
of Schirotzek et al. gives —0.64(7) fi+ |194j . The difference is probably mainly 
caused by an overestimation of the screening of the interaction, as also explained 
in Section 4.3. 1| We come back to this point in the next section. In Fig. 22 '&) 
the flow fjf_(l) of the spin-down chemical potential is depicted as a function of 
I. 
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5.4- Homogeneous phase diagram 

Next, we determine the nonuniversal critical properties of the strongly in- 
teracting Fermi mixture at nonzero temperatures with the RG approach [193] . 
and in particular the loca tion of the tricritical point in the homogeneous phase 



of is given by 



diagram. We use Eq. (78) for the flow of Vq , while the expression for the flow 



dj^ = A 2 /_ g (A) 

dA W-V^+Z'^A* ' [ ' 



similar to Eq. (85). The flow equation for Z' 1 can be obtained analytically 



from Eqs. (82) and (83), but is somewhat cumbersome to write down explicitly. 



Note that now we do not use the substitution ihuj m — > e q — /i+ anymore in 



Eq. (82), since this substitution is only correct for the extremely imbalanccd 
case at zero temperature. The initial conditions are the same as before with 
in addition fJ.' + (0) = As mentioned previously, the critical condition is the 
Thouless crite rion that the fully renormalized vertex V ' (oo) flows to zero. 
From Eq. (87), we see that incorporating the coupling constant Z' ~ solves the 



problem of the diverging chemical potential. 

We first apply the above procedure to study the equal density case (i.e., 
h 0) as a function of the negative scattering length a. The scattering length 
enters the calculation through the initial condition of V ' . To be able to 
express our results in terms of the Fermi energy ep = ep a , we calculate the 
(renormalized) densities of the atoms from 

dk 1 



(2tt) 3 e£( £ k-^(°°)) + 1' ^ 
In the weak-coupling limit, a — > _ , the chemical potentials hardly renormal- 



ize, so that only Eq. (78) is relevant. The critical temperature becomes expo- 



nentially small, which allows us to integrate Eq. (781 exactly with the result 
fceT c = 8epe 7 ~ 3 exp{— 7r/2/cF|a|}/7r and 7 Euler's constant. Compared to the 
standard BCS-result we have an extra factor of 1/e, coming from the screening 
effect of the bubble diagram that is not present in BCS theory. It is to be 
compared with the so-called Gor'kov correction |212j . which is known to reduce 
the critical temperature by a factor of 2.2 in the weak-coupling BCS-limit |213j . 
The difference with our present result is that we have only allowed for a nonzero 
center-of-mass momentum, whereas to get precisely the Gor'kov correction we 
would also need to include the relative momentum. We see that due to our 
approximation of neglecting the relative momenta in the bubble diagram, we 
overestimate the screening effect on the critical temperature by 20%. Note that 
in the previous section, we also already concluded that screening effects are 
slightly overestimated by the present simplest nontrivial RG. At larger values of 
|a|, the flow of the chemical potential becomes important and we obtain higher 
critical temperatures. In the unitarity limit, when a diverges, we obtain that 
T c = 0.14T F and fi(T c ) = 0.59e F . 
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Figure 23: The phase diagram of the homogeneous two-component Fermi mixture in the 
unitarity limit, containing the superfluid Sarma (S) and BCS phases, the normal phase (N) 
and a forbidden region (FR). The solid black line is the result of the RG calculations. The 
dots with error bars are experimental data along the phase boundaries as determined by Shin 
et al. [72| . The dashed and dashed-dotted lines are only guides to the eye. 



We can also calculate the critical temperature as a function of polarization 
p with the RG approach [193] and compare with the recent experiment of Shin 
et al. [75] . The result is shown in Fig. 23 The location of the tricritical 
point (TP) is determined from the fourth-order coefficient in the Landau theory 
for the superfluid phase transition. If this coefficient changes sign, then the 
nature of the phase transition changes from second order to first order. This 



was explained more elaborately in Section 3.4 where the right-hand side of Eq. 



(34) is the mathematical expression for the Feynman diagram in Fig. 24 'b), 
which is the diagrammatic representation of this fourth-order coefficient. In 
our RG approach, we calculate this one-loop diagram, where during the flow 
the integration over internal momenta is confined to the infinitesimal shell that 
is integrated out. Moreover, we have that during the RG flow the chemical 
potentials are renormalizing, so that self-energy corrections to the Feynman 



diagram of Fig. 24 b) are automatically taken into account. With this procedure 
we obtain for the tricritical point that P C 3 = 0.18 and T c3 = 0.063 Tp + . 

The results of Fig. [23] with the RG approach should be compared with those 



of Fig. 12 with the mean-field approach. We see that the location of the tricriti- 
cal point obtained from RG theory has shifted to much lower temperatures and 
polarizations than the one from mean- field theory. As a result, the RG calcu- 
lation is in much better agreement with experiments. We believe that the RG 
captures two main shortcomings of the mean-field theory, namely it takes into 
account fermionic self-energy effects and screening effects. Actually, the level 
of agreement with experiment is rather remarkable considering the simplicity of 
our RG. To some extent this is a coincidence, since there are many couplings 
whose renormalization we have ignored here although they could have a quan- 
titative influence, such as e.g. the effective mass of the fermions. In Ref. [193] . 
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we for example also included the center-of-mass frequency dependence of the 
interaction and found P C 3 = 0.24 and T c3 = 0.063 Tp + . Moreover, the results of 
the RG are also sensitive to the precise way in which we flow, so that the results 
depend for example on the intermediate cut-off HA' Q . We pick A' such that the 
high-energy two-body physics has been integrated out to a large extent, but the 
many-body physics not yet. This means that we take Aq to be a few times the 
Fermi wavevector. However, this procedure has some arbitrariness, and in an 
exact treatment the results should be fully independent of A' . We note that 
in Fig. [23j the dashed and dashed-dotted lines have the same meaning as in 
the homogeneous phase diagram of Fig. [l2j However, with our current RG for 
the normal phase these lines cannot be calculated, since for this a treatment of 
the supcrfluid phase would be required. Finally, we mention that at zero tem- 
perature, the Monte-Carlo treatment of Lobo et al. predicts a quantum phase 
transition from the equal-density superfluid to the polarized normal phase at 



a critical imbalance of p = 0.38, as was discussed in Section 4.3.2 |150j . This 
value seems to be in reasonably good agreement with experiments as seen from 
Fig. [23] 

5. 5. Renormalization theory for Cooper pairs 

In the previous paragraphs we studied renormalization effects using the 



fermionic action of Eq. (40). In Section 4.1 we showed that the Hubbard 



Stratonovich transformation can be used to derive two more actions that are 



equivalent to Eq. ( 40 1 . These two can also be used for a renormalization-group 



study of the strongly-interacting Fermi gas. In Refs. [208, 209] EH] EE] , the 



Bose-Fermi action of Eq. (47) was used as a starting point, while next we study 



renormalization effects for the Cooper-pair action of Eq. ( 52 ) . The Cooper-pair 
action is an exact microscopically derived action for the order parameter to our 
interest, which is the expectation value of the BCS pairing field A(r, r). In Sec- 
tion [3] we extensively studied the saddle-point or mean-field approximation to 
the Cooper-pair action. This simple approximation was seen to give a very satis- 
factory qualitative description of the experimentally observed phase transitions 
in the strongly-interacting regime. In Section [4~2] we also studied the Gaussian 
fluctuations of the pairing field in the normal phase for the balanced Fermi gas, 
called the NSR theory. This theory resulted in good quantitative agreement 
with thermodynamic measurements on the equation of state. Moreover, it was 
shown in Ref. [185] that the agreement between the Gaussian fluctuation theory 
and experiments is also good in the superfluid phase. 

Unfortunately, as soon as we consider even the smallest population imbal- 
ances, the Gaussian fluctuation theory gives rise to unphysical results [146] . 
Namely, near the critical temperature the Nozieres-Schmitt-Rink theory pre- 
dicts a negative polarization p = (n + — + n_) for a positive chemical 
potential difference — which corresponds to a compressibility matrix 
— d 2 Lu gc /d[i tT d[jL a -> with w gc the grand-potential density that is not positive def- 
inite. This is a very unsatisfactory situation, especially considering the success 
of the NSR theory for the balanced case. In this section we try to solve this fun- 
damental problem of the Gaussian fluctuation theory with the renormalization- 
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group approach. The NSR theory for the balanced Fermi gas [129J takes into 
account only a noninteracting gas of non-condensed Cooper pairs. We can 
therefore improve the theory for the Cooper-pairs by taking also into account 



the effect of the interactions between the pairs. Namely, as seen from Eq. (52), 
the pair action does not only contain a noninteracting part, but also two-pair 
interactions, three-pair interactions, and all higher-order interactions. A Popov 
theory for the bosonic pairs that includes the pair interaction effects was for- 
mulated by Pieri and Strinati, leading in the BEC regime to Popov's results for 
point-like bosons 216, 217] . In this regime the theory of Pieri and Strinati gives 
rise to a scattering length for the pair interaction of 0.7a with a the fermionic 
scattering length. This is rather close to the later obtained exact result of 0.6a 
|218j . Below the critical temperature, also the Bogoliubov theory for interact- 
ing Cooper pairs was studied [1221 11961 1183J . Other strong-coupling approaches 
that go beyond the NSR theory include so-called self-consistent ladder approx- 
imations [IBB1H7S] and Monte-Carlo calculations [2151 IW1 ITTH 11751 1220] . 

In this section, we apply the renormalization-group techniques developed in 
the previous section to study non-condensed interacting Cooper pairs in the 
unitarity limit. To this end, we have to generalize the Wilsonian RG theory 
for point-like bosons to the more complicated case of Cooper pairs |221j . The 
inverse propagator for the non-condensed Cooper pairs G^ 1 that follows from 



the quadratic part in the pairing field of Eq. (52 1 is given by Eq. (55 1. Here, 
we note that G~^ {ioj n ,k) only depends on the magnitude of k. The Feynman 
diagram that corresponds to the Cooper-pair propagator is shown in Fig. |24[a). 
We call G^ 1 the bare or microscopic propagator, indicating that no Cooper- 
pair interaction effects have been taken into account yet. Note that the bare 
propagator is exact, in the sense that it follows from an exact transformation of 
the fermionic action. With the RG approach we can consequently systematically 
include Cooper-pair interaction effects that lead to self-energy corrections to the 
bare propagator. 

The Cooper-pair interaction Va follows from the quartic part in the pairing 



field of Eq. (52) and is diagrammatically represented in Fig. 24 b) [216, 217 . 
Here, we do not take the full frequency and momentum dependence of the 
Cooper-pair interaction vertex into account, but we consider only two external 
frequencies and momenta to be nonzero, namely either u>i = — LU2 and ki = — k2, 
or w 3 = — w 4 and k 3 = — k 4 , where the labeling is given in Fig.[24p. This specific 
choice corresponds physically to considering only zero center-of-mass frequencies 
and momenta, which is motivated later on. The resulting expression is given by 

= wpv £ Go - (n '' k ') G °-+(- n '< - k ') 

n' ,k' 

xGo._(>' + n, k' + k)G 0)+ (-n' - n, -k' - k), (89) 

where we defined Va = Va(0, 0), so that Gy encapsulates the considered (rel- 
ative) momentum and frequency dependence of the Cooper-pair interaction. 
The function Gy is thus equal to the second and third line of Eq. [89] divided 
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Figure 24: Diagrammatic representation of a) the bare Cooper-pair propagator and b) the 
bare Cooper-pair interaction. The Cooper pairs are represented by thick lines, while the thin 
lines correspond to fermionic propagators. 



by Va- The Matsubara sum over odd fermionic frequencies n' in Eq. (89) can 
be performed analytically. Performing the sum for zero external momentum 
and frequency, we find for example the expression of Eq. (34). Indeed, the 
fourth-order coefficient /3l in the Landau theory can physically be interpreted 
as an interaction strength for the Cooper pairs. We call V&(iui n , k) the bare or 
microscopic interaction, in order to make the distinction with the effective or 
renormalized Cooper-pair interaction V A , which includes the effect of Cooper- 
pair fluctuations and that is calculated next during the RG flow. 

Due to the repulsive interaction between the Cooper pairs, they acquire a 
self-energy T, A . In the simplest approximation, the self-energy is momcntum- 
and frequency independent, so that the full renormalized propagator becomes 
G' A (iuj n , k) 1 = G A X (icon, k) — £a- We define the bare Cooper-pair chemical 
potential as /ia = ^G^ 1 (0, 0), while the renormalized chemical potential is given 
by fj,' A = /ja — ^Sa- The renormalized chemical potential thus includes pair 
self-energy effects. As a result, the full Cooper-pair propagator G' A (ioj n , k, /i A ) 
depends on the renormalized chemical potential /i A . The simplest RG calcula- 
tion that gives nontrivial results treats both the renormalization of the chemical 
potential fj,' A and the interaction strength V A . It ignores the three-pair interac- 
tions and higher. The flow equations for these two coupling constants can be 
derived along exactly the same lines as for point-like fermions, as we did in Sec- 



tion 5.1 The corresponding one- loop Feynman diagrams are diagrammatically 



represented in Fig. 25 Note that a thick line denotes a Cooper pair propagator, 
whereas the thin lines in Fig. [24] denote fermion propagators. Since we are now 
performing a RG study solely for the Cooper pairs, only the momenta of the 
thick lines are restricted to the considered high-momentum shell. The bare cou- 
pling constants for the pair action are derived from the Hubbard-Stratonovich 
transformation and integrating over all fermionic fields. As a result, the inte- 
gration over momenta in Fig. 24 for the thin lines is thus in the present case 
over all momenta. 

When doing the derivation of Section |5~l] also for the Cooper pair RG equa- 
tions, a few things should be kept in mind. The Cooper pairs do not carry 
spin as the fermions, so that the Cooper pair interaction is between equal pairs, 
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Figure 25: Diagrammatic representation of the '/3 functions', a) Feynman diagram determing 
the self-energy of the Cooper pairs, b) Feynman diagrams renormalizing the Cooper-pair 
interaction. The middle diagram is also called the ladder diagram, the right diagram is 
called the bubble diagram. Note that the lines are thick and correspond to the Cooper-pair 
propagator. 



in contrast to the interaction for fermions between opposite spin. As a result, 
the Cooper pair interaction parameter is divided by the usual factor of two in 
the pair action to avoid double counting problems. Another difference with the 
fermionic action is that the frequency dependence of the Cooper-pair propagator 
is more complicated. As a result, the Matsubara sums of the one-loop Feynman 
diagrams cannot be performed analytically anymore, but have to be evaluated 
numerically within each momentum shell. More detailed accounts for deriving 
the RG equations for bosons can be found in Refs. [1271 I5U] . 

The RG flow of the Cooper-pair interaction strength and chemical potential 
are found to be determined by the following set of coupled differential equations 



dfA 
dl 



dVk 
dl 



where the '/3-functions' are given by 



ft 



A, 



n 

2vr 2 



{H(A,,/4)+4II(A,,/4)}. 



(90) 



(91) 



(92) 



Here, II and S are the so-called 'bubble' and 'ladder' contributions to the effec- 
tive Cooper-pair interaction, which we also encountered in the RG theory for 
the fermionic action. Moreover, A; denotes the wavevector of the Cooper pairs 
in the shell of infinitesimal width. This wavevector is parametrized by the flow 
parameter I, and we start the RG flow at the high- momentum cutoff HAq and 
decrease as A; = Aoe~'. In addition, A/ is the derivative of A; with respect to I. 



Solving Eq. ( 90 1 for increasing / means that we are including the effect of pair 
fluctuations with lower and lower momenta, while due to the coupling of the 
differential equations we automatically generate an infinite number of Feynman 
diagram s, s howing the nonperturbative nature of the RG. The initial conditions 
for Eq. Q are fi' A {l = 0) = ha = ^G A 1 (0,0) and V A (l = 0) = V A , which are 
calculated from Eqs. (55) and (89 1. 
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The one-loop expression for the renormalization of the chemical potential in 



Eq. (91 1 that determines the self-energy of the Cooper pairs, has a clear physical 
meaning, since it is seen to be proportional to the renormalized pair interaction 
strength and to the density of Cooper pairs. The 'bubble' diagram II(fc, /i' A ) 
describes the effect of 'particle-hole' excitations on the effective Cooper-pair 
interaction, where these 'particles' are now actually Cooper pairs. It is given by 

n(k, t i' A ) = jJ2G'A(iun,k,ti' A ) 2 . (93) 

n 

The 'ladder diagram' describes the Bose-enhanced scattering of the bosonic 
Cooper pairs, given by 

E(*,Ma) = jJ2 G y( iuj n,kf\G' A (iuj n ,k, f ,' A )\ 2 , (94) 



where the momentum and frequency dependence of the interaction in Eq. ( 89 ) 



Gy, is seen to enter. This frequency and momentum dependence is im- 



portant, since otherwise Eq. (94) ultimately would lead to an ultraviolet diver- 
gence. The divergence physically arises from approximating the pair interac- 
tion as a point interaction, which is therefore insufficient. We also note that 



the self-energy diagram of Fig. |25[a) and Eq. (91), and the bubble diagram of 



Fig. 25 'b) and Eq. (93) do not lead to divergencies. As a result, our present 
scheme for including the Cooper-pair interactions is the minimal choice for ob- 
taining divergence-free, or equivalently, cutoff independent results. 

5. 6. Results 

In Section |4.2| it was explained that the NSR theory of the strongly inter- 
acting normal state gives rise to two contributions to the grand-canonical ther- 
modynamic potential, namely a contribution describing an ideal gas of fermions 
and a contribution describing an ideal gas of noncondensed Cooper pairs |129j . 
The contribution to the grand potential density due to the Cooper pairs is given 
in our RG approach by the one-loop expression 

=- A <2^E lo g[- G A~^W A )L (95) 

n 

where the first minus sign on the right-hand side is only present when A/ is a 
decreasing function. Note that this last expression reduces to the differential 
form of the NSR contribution to the grand potential density in Eq. |58| ) when 
the Cooper-pair chemical potential is not renormalized (pf A = (J.a), i-e., when 
we consider the Cooper pairs to be noninteracting. If the exact Cooper-pair 



propagator is inserted in Eq. ( 95 ) , then the exact grand potential density is 
obtained. However, this would require the treatment of all n-body interactions, 
which is presently out of reach. 
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To be able to evaluate Eqs. (90) and (95) numerically, it is convenient to 



perform contour integration. The Green's function of the Cooper pairs from 



Eq. (55) can be written in the spectral form [222 



G A («w„,fc,^ A ) = - / duj ^ ; 

7T ./ W — lUJ n 



(96) 



where lo + = uj + ir\ with r\ \ 0. The imaginary part of the Green's function 



can be obtained analytically and was given in Eq. (59). With the spectral 



representation, we can rewrite Matsubara sums over the pair Green's function as 
frequency integrals that are convenient for numerical evaluation. For example, 
we have 



^E G A( lWn ' fc '' 1 A) = ^ J dujb ^) lm \ G '&{u + ,k,fj,' A )] 



(97) 



where n is even and 6(e) = l/(e^ e — 1) is the bosonic distribution function. 

(98) 



Moreover, the pair bubble diagram from Eq. ( 93 ) becomes 

n(fc, ju A ) = ^2 G' A {iuj n , k, Ma 
2 



hp 

cLob(u)Im[G A (<j + ,k,ti' A )]Re[G' A (u; + ,k,Li A )}, 



where we used Eq. ( 96 ) and the Kramers-Kronig relation to relate the real and 



imaginary part of the Cooper-pair Green's function. For the grand potential 



density from Eq. ( 95 ) , the frequency integral form is given in Eq. ( 58 ) 



After having performed the RG calculations, we have that the total thermo- 
dynamic potential density is given by 



(T,fi a ) = ^ T ^^ = Wlg (T,/v) +WA,oo(r,M CT ), 



(99) 



with ft the grand potential, V the volume, and cja.oo (T, fi a ) = ^a(1 — > oo). 
From Wgc, all other thermodynamic quantities of interest can be obtained by 
the standard thermodynamic relations. For the balanced case, the results of 
the renormalization procedure are shown in Fig. |16| Here, it is seen that due 
to the inclusion of the repulsive Cooper-pair interactions the equation of state 
shifts slightly down towards the noninteracting result. However, the effect is 
not strong and the agreement with the experiment is still reasonable. We note 
that a more elaborate treatment of fermionic self-energy effects on the Cooper 
pair propagator would increase the pressure again [175| 1185) . The results of 
the present procedure for the imbalanced case are shown in Fig. |16| Here, we 
have calculated the equation of state for the imbalanced Fermi gas. In Fig. [26j 
we show the pressure p g as a function of h/fJ, for the temperature T '/ ' \i = 0.75, 
where h = (fi + — /i_)/2. At this temperature, the NSR approximation predicts 
a negative polarization for positive h, which is an unphysical result. The present 
RG approach that treats the Cooper pair interaction effects beyond NSR theory 
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Figure 26: Equation of state for the normal phase of a strongly interacting imbalanced Fermi 
gas at unitarity. The pressure p g = — ui gc (T, fx, h) of the gas is calculated at temperature 
T = 0.75/1 as a function of h//i = — fjt—)/([j,+ + fi—) with the renormalization group 
approach (full line), the Nozieres-Schmitt-Rink approach (dashed line), and for the ideal Fermi 
gas (dashed-dottcd line). For each curve, the pressure of the imbalanced gas is normalized 
to the corresponding pressure of the balanced gas po = — 6j gc (T, fi, 0). We note that po is 
different for the three different methods, as can be seen from Fig. |16| 



does not have this problem. We see that, as a result, the present RG theory, 
the NSR theory and the mean-field theory give very different results for the 
pressure of the imbalanced Fermi gas in Fig. [26| This pressure has recently 
been measured at zero temperature, where good agreement with Monte-Carlo 
calculations was obtained [75] , 

Finally, we discuss the effect of the Cooper-pair interactions on the critical 
temperature for the balanced Fermi gas. In the unitarity regime, the effective 
two-pair interaction is repulsive in the normal state. Taking into account the 
repulsive two-pair interaction makes the ratios k-QT c / fi and £f/m smaller, be- 
cause the repulsion lowers the effective chemical potential of the non-condensed 
Cooper pairs. As a result, the density of non-condensed Cooper pairs becomes 
lower, which decreases the total density and the critical temperature for con- 
densation. More precisely, when G^ 1 (0,0) = /ia = 0, then both mean-field 
theory and NSR theory predict a transition to the superfluid state, i.e., the con- 
densation of Cooper pairs at k B T^ { = 0.66/x. Upon lowering the temperature 
from T^£, jUa becomes positive, but due to the repulsive interactions the renor- 
malized chemical potential goes down, and at the end of the flow /Lt' A ^ might 
still be negative. Therefore, the critical condition in the presence of Cooper- 
pair interactions becomes ^ = 0. With our present approach, this results 
in fce^rg — 0.43/4, which is much closer to the most recently measured value 
of fceT^p = 0.40/j, [173] . Note that this result could be further improved by 
performing a RG calculation for the superfluid state. Namely, close to T c we 
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then initially have fJ.' A {l) > 0, and the RG flow starts out in a superfluid state. 
Due to the repulsive interactions, the effective chemical potential is then again 
lowered, and the system can flow into the normal state |127j . 

Having discussed in this section in more detail the effect of interactions on 
strongly-interacting Fermi mixtures for the homogeneous case, we look next 
more accurately at the effects of the trapping potential, which we have so far 
only treated in the local density approximation. 



6. Inhomogeneous Fermi gases 

Ultracold atom experiments arc always performed in a trap to avoid contact 
of the atoms with material walls that would heat up the cloud. Due to this 
trapping potential the atomic cloud is never truly homogeneous. However, typ- 
ically the trapping frequency corresponds to a small energy scale in the system, 
so that the inhomogeneity is not very prominent. In this case, we may use the 
so-called local-density approximation (LDA), which was discussed in Section 



3.5| It physically implies that the gas is considered to be locally homogeneous 
everywhere in the trap. The density profile of the gas is then fully determined 
by the condition of chemical equilibrium, which causes the shape of the cloud 
to follow an equipotential surface of the trap. 

But even if the trap frequency is small, LDA may still break down in a few 
specific cases. A first example is near the edge of the cloud, where the local- 
density approximation predicts at zero temperature that the particle density 
goes strictly to zero beyond the radius where the chemical potential equals 
the potential energy in the trap. This is not fully correct, because quantum- 
mechanically there is an exponentially suppressed probability that the particles 
reside in the (semi-)classically forbidden areas. Another example occurs when 
an interface is present in the trap due to a first-order phase transition. For a 
resonantly interacting Fermi mixture with a population imbalance in its two 
spin states j24j (25] , such interfaces were encountered in the experiments by 
Partridge et al. [35] and by Shin et al. [75] at sufficiently low temperatures. 
Here, the application of LDA leads to a discontinuity in the density profiles of 
the two spin states, which would cost an infinite amount of energy when gradient 
terms are taken into account. Experimental profiles are therefore never truly 
discontinuous, but are always smeared out. 

The presence of a superfiuid-normal interface also can have further conse- 
quences. Namely, in a very elongated trap Partridge et al. observed a strong 
deformation of the minority cloud at their lowest temperatures. At higher tem- 
peratures the shape of the atomic clouds still followed the equipotential surfaces 
of the trap [71] . A possible interpretation of these results is that for tempera- 
tures sufficiently far below the tricritical point I Hi. 85], the gas shows a phase 
separation between a balanced superfluid in the center of the trap and a fully 
polarized normal shell around this core. Because the trap shape is ellipsoidal, 
the superfluid core could consequently be deformed from this shape due to the 
surface tension of the interface between the two phases, which prefers shapes 
with the least possible area [149 . This deformation then causes an even more 
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dramatic break down of LDA. We note here that to understand the large surface 
tension of the Rice experiment that causes rather extreme deformations in their 
case, it is probably needed to take non-equilibrium effects into account |151| [78] . 

In order to further investigate the inhomogeneity of trapped Fermi gases, we 
thus need to go beyond the local density approximation. There are several ap- 
proaches to incorporate the inhomogeneity of the system beyond LDA, of which 
we discuss two frequently used methods. First, we discuss an approach based 
on Landau-Ginzburg (LG) theory, which includes a gradient energy penalty for 
a spatially varying order parameter. We try to incorporate the knowledge from 
normal and superfluid Monte-Carlo equations-of-state into the LG functional 
and calculate the surface tension of the superfiuid-normal interface in the uni- 
tarity limit. We also compare the LG approach with the interface observed at 
the MIT experiment [72]. Next, we discuss the Bogoliubov-de Gennes (BdG) 
equations, which take the single-particle states of the trapping potential into 
account. As a result, the BdG approach gives the exact result for the noninter- 
acting trapped Fermi gas, so that also the edge of the noninteracting cloud is 
correctly described. The interactions between the particles in the BdG approach 
are treated on the mean-field level, so that the results from these equations can 
be trusted at best qualitatively. We will see that the BdG equations give rise 
to oscillations in the superfluid order parameter and the particle densities near 
the superfiuid-normal interface. These oscillations are interpreted in terms of 
the so-called proximity effect. 



6.1. Landau-Ginzburg approach 

When using Landau-Ginzburg theory, we write the grand-canonical ther- 
modynamic potential as an expansion in the order parameter, where usually 
only the first few local terms and gradient terms are kept, as was for example 
done in Eq. (37). In most cases, it is not possible to find exact expressions 
for the coefficients in the expansion, so either mean-field theory, or other ap- 
proximations have to be invoked in order to determine these parameters. As we 
discussed in Section [3] the BCS mean-field theory seems to give a satisfactory 
description of the qualitative phase diagram for the imbalanced Fermi gas, as 
determined by recent experiments [72] • We therefore use the BCS theory as a 
starting point for our discussion of the inhomogeneous Fermi gas. In Section 
|3.6[ we discussed the phase diagram for the trapped Fermi gas, where we used 
the local-density approximation to account for the trapping potential. This 
means that we applied the homogeneous theory with spatially varying chemi- 
cal potentials, /v(r) = \i a — V cx (r) = [i a — mu) 2 r 2 /2. Due to the fact that a 
first-order transition could occur as a function of the chemical potentials, the 
LDA approximation predicted a jump in the superfluid order parameter and the 
corresponding particle densities. 

As just mentioned, a true jump in the order parameter is unphysical, and the 
easiest way to go beyond LDA for a more realistic description of the superfiuid- 
normal interface is to take also into account the gradient energy of the order 
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parameter in the grand potential functional by 



fi[A;M] = / dr { 7 |VA(r)| 2 + ^ gc [A(r); (x(v), h}} . (100) 

with w gc [A;/x, h] the local grand potential density. We encountered the coeffi- 
cient 7 already in Section |3.7[ where 7 was interpreted as being inversely pro- 



portional to the effective Cooper pair mass. Note that in Eq. ( 100 ) we neglect 
the dependence of 7 on A and we also ignore all terms that are of higher order 
in the gradients, resulting in the simplest extension beyond the local-density 



approximation. In Section 3.7 we showed how to calculate 7 from the normal 
state using the Cooper-pair propagator, diagrammatically represented in Fig. 
[24} A second way to compute a value for 7 is to use the fact that this coefficient 
can be related to the superfiuid density p B . Namely, it is a standard result from 
Landau-Ginzburg theory that 7 = h 2 p s /8m 2 \ (A) | 2 |223j . Using the knowledge 
from the superfiuid Monte-Carlo equation of state, we get the zero-temperature 
result 



where we used that at zero temperature the superfiuid density is equal to the 
total density. We have that /3 MC = -0.58 [2H)I [Wl |2"2"4] . which was already 
encountered in Section |4.3.2| and that Cmc = 1-07, which gives the Monte- 
Carlo value [225] for the universal relation between the pairing gap and the 
chemical potential for the unitary Fermi gas, (A) = Note that by using 
BCS theory, we would have obtained 7bcs = V^/^^^Cbcs^ + Pbcs) 3 ^ 2 V^^ 
where /3bcs = —0.41 and Cbcs = 1-16 are the mean-field results for the same 
universal constants. 

6.1.1. Surface tension 

Since we are studying the superfluid-normal interface beyond the LDA, we 
can determine the surface tension of the calculated interface |144j . The surface 
tension is the work per unit area that has to be done to increase the area of the 
interface. Since the LDA assigns no energy to the interface, we can calculate 
the surface tension from the difference in the grand potential between the LDA 
result with a discontinuous step in (A)(r) and the Landau-Ginzburg result with 
a smooth profile for the order parameter (A) (r) . We will calculate the surface 
tension by considering a flat interface in a homogeneous box, which has an 
interface around z = with a normal phase at the bottom of the box (z < 0) 
and a superfiuid phase on top (z > 0). We do not consider any gravitational 
effects on the gas. The superfluid-normal interface occurs when the imbalance is 



critical, i.e., when h — h c (fi) = n/i, where in Section 3.4 we found kbcs = 0.81 
using mean-field theory, while in Section |4.3.2| we found kmc = 0.93 using 
the results from the Monte-Carlo equation of state |150) . Since at the critical 
imbalance the grand potential of the normal state minimum is exactly equal to 
the grand potential at the superfiuid state minimum, we have that the surface 
tension is given by the difference between a system that stays in one minimum 
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and a system that goes near the interface smoothly from one minimum to the 
other. 

To determine the resulting shape of the interface in the box, we use the 
grand potential 

n[A;(i,h c ] =Ajdz j 7 ^Mj +u> sc [A(zy,v,h c ]X. (102) 

where we condider a real-valued order-parameter profile A(z), and A is the area 
of the box at z — 0. Extremizing this functional results in the Euler-Lagrange 
equation for the interface 

du} gc [A;^h c ] d 2 A 

M 27 d^ = - (103) 

which can be readily solved numerically. This results in a smooth monotonic 
function (A) (z) that on the normal side of the interface approaches zero and on 
the superfluid side approaches the equilibrium value of the superfhiid minimum 
(A). Inserting this profile in Eq. (100) and computing the difference with 
S7[0; fi, h c ] then results after division by the area in the surface tension. 

There is also a more elegant way to compute the surface tension a, which 
does not explicitly require the shape of the interface 226 . Namely, we have 
that 

poo 

dz (u} m [(A)(z);n,h c ] -u} m [0;n,h c }), (104) 

which can be rewritten as an integral over a new integration variable A', using 
that (A)(z) is a monotonically increasing function between zero and (A). Doing 
so, we obtain 

r {A) dz 

° = J d ^'dA' ( w sc[ A ';mA] -w gc [0;MAD, (105) 

which still requires the determination of dz/dA' . This can be done by consid- 
ering 

d fdA\ 2 dAd 2 A ldAduj RC , . 

' - = = -^T^O (106) 



dz \ dz J dz dz 2 7 dz dA ' 



where we used Eq. (103). After integration over z this gives 
dA\ 2 1 



dz j = - (w g c[A(z);/i,ft, c ] -(J sc [0;^,h c ]), (107) 

where we used that (dA/dz) 2 = at z = —00. Upon taking the square root of 
Eq. (107), inverting the result and insertion in Eq. (105), we obtain the final 
result for the surface tension 

"(A) 



a 







dA'^7 (w gc [A' ; /1, h c ] — ^gc [0;»,h c ]), (108) 
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which does not explicitly depend on the solution (A) (z) anymore, but only on 
7 and the shape of the barrier in the grand potential density uj gc . Writing the 
surface tension as a = rj^m/H 2 , we have that the dimensionless number r] only 
depends on the temperature. So in order to calculate the surface tension we 
need an expression for the grand potential as a function of the order parameter 



A. Using the BCS mean-field thermodynamic potential from Eq. (27), we find 
that 7/bcs = 8.6 x 1CT 3 . 

6.1.2. Monte- Carlo improved grand potential 

In Section |4j we discussed the results of Monte-Carlo calculations for the 
superfluid equation of state and the normal equation of state in the unitarity 
limit. There we saw that the BCS mean-field theory gives in particular a rather 
poor description of the strongly-interacting normal state. In Section 4.3.1 we 
also proposed a simple way in which the Monte-Carlo equation of state for 
the normal state could be included in the grand potential by using the Ansatz 



for the 'renormalized' chemical potentials in Eq. (67). Here, we apply a sim- 
ilar procedure to the superfluid state. The first step is to insert the normal 
state renormalized chemical potentials in the BCS grand potential density, giv- 
ing Wi mp (A;/v) = w bcs(A; so that the strongly-interacting normal 
equation-of-state is automatically included when A = 0. 

Although the Monte-Carlo results for the normal state are now included, 
the grand potential needs to be further modified in order to also give agreement 
with the Monte-Carlo superfluid equation of state |224] . for which we have 
at zero temperature that fi = 0.42eF and (A) = 0.45ef = 1-07^ [225 . One 
issue with the construction so far is that <9u>i mp ((A); /i, h)/dh is nonzero in 
the superfluid state, because (i9wbcs((A); //, h')/d/j,')(dij' /dh) is nonzero. The 
reason for this behavior is that the self-energy corrections captured by // depend 
on h, which is correct for the normal state, but not for the superfluid state, 



where the experiments of Section 4.3.2 show an equal-density superfluid core in 



the unitarity limit, so that the grand potential should become independent of 
h. This behavior is correctly captured by the BCS potential itself without the 
'renormalization' of the chemical potentials, which indeed becomes independent 
of h for A > h. From this discussion, we conclude that // should become 
independent of h as a function of A. Another issue is that the parameter A in 



Eq. (67), which describes the interaction effects in the normal state, is too big for 
the superfluid state, since in the latter case the BCS mean-field theory already 
takes into account a large part of the interaction effects through the presence 
of the Cooper-pair condensate A. We thus expect an improved result if we let 
the self-energy parameter A decrease as a function of |A| 2 /^ 2 . Assuming that 
there is no h dependence of \j! in the superfluid state, we find that by using 
A' = A - B\A\ 2 /[i 2 we simultaneously satisfy fx = 0.41eF and (A) = 1.06/i, if 
we take B = 0.135, while A = 1.01 is the earlier obtained value for the normal 
state. 

However, to have a complete grand potential we still must specify how fjf 
evolves as a function of A in order to lose its dependence on ft, for A > h' . 
To ensure a smooth interpolation between the normal state and the superfluid 
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Figure 27: Upper panel: thermodynamic potential density ui for the unitary Fermi gas as a 
function of the superfiuid order parameter A. The dashed line gives the result for BCS mean- 
field theory, while the full line effectively includes the superfiuid and normal equations of state 
from Monte-Carlo calculations. Lower panel: the 'Monte-Carlo' improved thermodynamic 
potential density for the balanced case, h = (full line), for the critical chemical potential 
difference, h c = 0.99(j, (dashed line), and for h = 1.2h c (dash-dotted line). 



state, we use the switching function 
/.(A, ft') 









a- 


- ^cos 











d(h'-\A\) 



(109) 



which is a function of |A| 2 that equals 1 for A = and which goes to zero with 
vanishing derivative at |A| = h', the point where the modified BCS thermody- 
namic potential loses its dependence on ft'. The switching function stays zero at 
larger values of |A| due to the presence of the stepfunction 8. If we include the 
switching function in the expression for p! as a function of /i and ft, obtained 



from Eq. ( 67 1 , we get 
•5// 



/i 



ti = 



10 - 3A' 

5ft 
5 + 3 A' ' 




3(10- 3A')A' /ft 
(5 + 3A') 2 U 



f s (A,h')\, (HO) 



(111) 
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where we have also replaced A with A'. The latter is a function of |A|//x in 



the way mentioned above. The expressions in Eqs. (110 1 and (111) are valid 
for /i > h(l — 3A/5)/(l + 3A/5), while for smaller values of \x the system 
is in a fully polarized phase described by an ideal gas of majority particles. 



When A = 0, we have A' — A and / s — 1, so that then Eqs. (110) and 



(111) contain exactly the same information as Eq. (67). Note that the chosen 
switching function is arbitrary in the sense that it is solely motivated by its 
smooth interpolation properties between and 1. However, by including it we 
have obtained a grand potential that not only has the Monte-Carlo equation-of- 
state for the normal state and the equilibrium supcrfluid state incorporated, but 
also smoothly interpolates between these two regimes. We will use this zero- 
temperature 'Monte-Carlo improved' grand potential density ajj mp to compare 
with the results from BCS mean-field theory and experiments. It is plotted in 
Fig. [27| where in the upper panel Wi mp and wbcs are compared. The 'Monte- 
Carlo improved' potential is seen to give rise to a much larger pressure of the 
gas p g (p g = — w gc ), which is due to the effective inclusion of the self-energy 
effects. In the lower panel, we show u>i mp for different values of h, showing that 
at zero temperature the critical chemical potential difference is determined by 
h c = 0.99/i, which is very close to the Monte-Carlo result of h c = 0.96/1 150 . 

6.1.3. Results 

We calculate the surface tension with the use of Eq. (108), where we insert 
for the grand potential density Wi mp (A; /i, h) — wbcs(A; //, h!) with cjbcs, A 4 ' 



and h! given by Eqs. (27), (110) and ), respectively. We find that at zero 



temperature 77i mp = 1.6 x 10 _:i . This result is to be compared with the value of 
r/Bcs = 8-6 x 10 -3 from mean-field theory, and r)£ t = 0.6, which was fitted to the 
large deformations of the superfluid core observed in the Rice experiment [149] . 
The experiment at MIT of Shin et al. does not show any deformation, putting an 



upper bound on ?/ of about 0.1 |149j . This is in agreement with the above value 
of T/imp at zero temperature. For higher temperatures the surface tension would 
further decrease due to a lowering of the superfluid barrier, and ultimately vanish 
at the tricritical point. Understanding the large value for the surface tension 
found in the early Rice experiments is not possible within the present approach, 
since in that case non-equilibrium effects seem to be important [73 1151] . The 
experimental results of Shin et al. were shown already in Fig. |20| where they 
were compared with the theoretical results from the Monte-Carlo equations-of- 
state in combination with the LDA. In Fig. [28j the same experimental data is 
shown again, but now compared with the present LG approach. For simplicity, 
we have assumed here a spherically symmetric trapping potential, although the 
actual experimental trapping potential was elongated in one direction. The 
trapping frequencies for the MIT experiments were given in Section |3.5| As 
explained in the same section, the elongation can be scaled away by a coordinate 
transformation as long as the different directions in the trap are not coupled, 
as was done for example in applying the local-density approximation. Then, 
the elongation plays no physical role and by transforming the coordinates a 
spherically symmetric potential is obtained. However, in the Landau-Ginzburg 
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r/R + 

Figure 28: The density profiles of a unitary Fermi mixture in a harmonic trap with a 
spin polarization of P = 0.44. To realize this polarization, the chemical potentials = 
81hui, fi— = I9h£j, and temperature k^T = lOhui were used with Hi the (effective) trapping 
frequency. The upper figure shows the majority (upper blue curve) and minority (lower red 
curve) densities as a function of the position in the trap. The lower figure shows the density 
difference. The horizontal axes are scaled by = y / 2/i+/™ 2 , while the vertical axis is 
scaled by the central spin-up density no = n_|_(0). The experimental data (dots) are from 
Shin et al. [72]. 
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functional of Eq. ( 100 1 the directions in the trap are coupled by the gradient 
term, so that the mentioned coordinate transformation then actually transfers 
the effect of the trap anisotropy to the gradient term. It was shown in Ref. 
|226| that the effect of this anisotropy is actually small, which also agrees with 
the experimental results of Shin et al. that find a gas cloud which follow the 
equipotential lines of the trap very closely. To get the order-parameter profile 
for (A)(r), we solve the Euler-Lagrange equation for the spherically symmetric 



case with 7 determined by Eq. ( 101 1, where we used for /1 the local value in the 
trap at the interface. 

To obtain the density profiles, we use the local thermodynamic relation 
ncr(r) = — duj gc (T, /i CT (r))/<9/i CT , by which we mean that we differentiate the 
used grand potential density w gc = Wi mp with respect to the chemical poten- 
tials [i a and evaluate the result at the local value for the chemical potential 
fi a (r) = \i a — V ex (r) with V cyi (r) = mcD 2 r 2 /2 and the order parameter (A)(r). 
This allows for a full comparison between theory and experiment as shown in 
Fig. [28] The upper panel shows the majority and minority densities as a function 
of the radial position in the trap, while the lower panel shows the density differ- 
ence. The agreement is seen to be very good, although the data is somewhat too 
noisy for an accurate study of the interface. Although not very clearly visible in 
Fig. [28] we find that there is a small kink in the majority density profile. This 
kink signals the transition to the gapless Sarma phase. At the transition, the 
order parameter (A) becomes locally smaller than the renormalized chemical 
potential difference h! and the unitarity limited attraction is no longer able to 
fully overcome the frustration induced by the imbalance. As a result the gas 
becomes a polarized gapless Sarma superfluid. We found in Section [3^4] that for 
the homogeneous mass-balanced case at zero temperature, the Sarma phase was 
unstable. However, because of the inhomogeneity of the trapped gas, the order 
parameter is forced to move from its equilibrium value at the center of the trap 
to zero at the edge of the cloud. Therefore, at the interface we unavoidably en- 
counter locally a stabilized Sarma state, even at zero temperature. Notice that 
this is a feature of the smooth behavior of the gap and thus does not depend 
on the quantitative details of the grand potential density functional. 

6.2. Bogoliubov-de Gennes approach 

In the previous section we presented the simplest theory beyond the local- 
density approximation, namely a Landau-Ginzburg theory that takes into ac- 
count the gradient energy of the order parameter. Although this resulted in 
a nonlocal treatment of the order-parameter profile, the underlying fermionic 
wavefunctions were still treated in a local-density approximation as follows from 
the use of the homogeneous thermodynamic potential density with a spatially 



varying chemical potential in Eq. (100). To overcome this local approximation 



for the fermions, we now use a different approach to study the inhomogeneous 
system, namely the Bogoliubov-de Gennes (BdG) equations |156j . The BdG ap- 
proach has been widely used in the superconductivity literature to study a wide 
range of inhomogeneous situations, for example involving interfaces between su- 
pcrfluids and metallic or magnetic phases. The method has also been applied by 



several groups to study the supcrfluid-normal interface in an imbalanced Fermi 
gas [223 1213 [113 El Ell [215]. Some of these studies found that near the inter- 
face the order parameter becomes oscillatory, which was suggested to be related 
to the Fulde-Ferrell or Larkin-Ovchinnikov (FF or LO) phases [223 QH H3] 
that were discussed in the introduction. In a follow-up study [230] the authors 
of Ref. |142) showed that the oscillations are not a finite size effect |229) . but 
rather caused by the so-called proximity effect [2311 1232L 12331 1152j . In this sec- 
tion, we briefly review the Bogoliubov-de Gennes method, where we start with 
a derivation, after which we present an efficient numerical scheme to solve the 
obtained coupled di fferential equations. We apply the method to the phase- 
separated imbalanced case and indeed find the mentioned oscillations near the 
interface. Our findings support the conclusion that they are caused by the 
proximity effect. 



6.2.1. Derivation 

To derive the Bogoliubov-de Gennes equations, we start with the Hubbard- 
Stratonovich transformed action of Eq. (47). The action can be written in a 



2x2 matrix form, which we call spin or Nambu space, as we saw explicitly in 
Section |4.1[ As a result, we have for the quadratic part in the fermionic fields 
of the action that 



S (2) [<M 



dr J dr (0;(r,r),0_(r,r)) 

(hd r -h + H Q {v) A(r) 

^ A*(r) hd T -h-H (r] 



(112) 



-(r,r) 
-(r,r) 



where we have made the approximation that the pairing field is independent of 
(imaginary) time. The diagonal part of the 2x2 matrix contains 



H (r) = -— V 2 + V ox (r)- M , 
2m 



(113) 



where the external potential V cx (r) is added to describe the trapping potential. 
We consider a harmonic trapping potential, which approximates the experimen- 
tal situation very well. Just like in the previous section, we use a spherically 
symmetric trapping potential V c *(r) — mui 2 r 2 /2 with Q the (effective) trap 
frequency. 

Next, we perform a Bogoliubov transformation with the goal to diagonalize 



the position-dependent part of the 2x2 matrix in Eq. (112). To this end, we 
apply the following unitary transformation [156] 



Mr,r) 



E 



«n(r) -<(r)\ / il>+, n (r) 
v a (r) <(r) J U* „(t) 



(114) 



where ±n = (n,£,±me) denotes the set of three quantum numbers required 
to specify the single-particle eigenstates. We see that the unitary transforma- 



tion indeed diagonalizes the spatial part of Eq. (112), if the time-independent 
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Bogoliubov-de Gennes equations are satisfied, which are given by 

H Q (r) A(r) \ ( u n (r)\ _ fu n (r)\ 
A*(r) -H (r)J {v n (r) J ~ ^ [v n (r) J ' ^ 

This is a set of coupled differential equations, for which we still need to specify 
the boundary conditions and the self-consistency relation for A(r). The coher- 
ence factors u n and v n are normalized for each n, namely J dr (|w n | 2 + |w n | 2 ) = 1. 



Notice that when (u n , v n ) is a solution to Eq. ( 115 1 for energy E n , then (— u*, u* 
is a solution with — E n . Since the considered trap is spherically symmetric, the 
gap A(r) is a function of the radius only. We therefore have that 

u»{r,e,<f>) = u ^Y lme {e,<t>), (H6) 

r 

where Yg me are the spherical harmonics, and we do the same for v. The sum 
over n is now a sum over the set {n,£, mi} with I > and me = —£, ...,£. In 
the following, we will consider without loss of generality that u n i(r), v n e(r), and 
A(r) are real-valued functions. As result, we find the following equation for the 
functions u„e(r) and v n e(r), namely 

~H rA r)M r) ) • (117) 
dr 2 \Vni[r) J \Vne\r) ) 

The matrix H n e is given by 

H M- 2m (»-Vr(r)+Eni -A(r) \ 

Uni{r) -W[ A(r) /i-V 4 (r)-£^J' (U8) 

with V g ex (r) the effective external potential, i.e., 

which includes the effect of the centrifugal barrier for nonzero angular momen- 
tum I. The boundary conditions for u n i> and v n i are that they should be both 
zero in the origin and at infinity. 

Before we start with the discussion of the numerical procedure to solve the 
BdG equations, we give the analytic solution for the non-interacting case, i.e., 
A(r) = 0. Then, the differential equations are not coupled, so that they reduce 
to the 3-dimcnsional harmonic oscillator problem, whose solution can be found in 
many textbooks. The properly normalized noninteracting solutions (w n , v n ) = 
(</Wm,0) are given by 

&am(r,M) = Mnte-^ L e n +l " (^pjY em (e,<j)), (120) 
Ent = tiw(^+2ri + £\ -it, (121) 



/ 2 n + e + 2 n ] 

Mnl = VP(l + 2£ + 2n)\\V^' (122) 
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with L l n the associated Laguerre polynomials, M n t the corresponding normaliza- 



tion constants, and I = ^Jh/mui the so-called harmonic trap length. Eq. (1151 
also gives rise to the solutions (u n ,v n ) — (0, 4> n im) with energy —E n g. 

6.2.2. Numerical methods 

There are multiple ways to solve the second-order matrix differential equa- 



tion of Eq. ( 117 1. For example, it is possible to use the (renormalized) Numerov 
method, which is a well-established numerical approach for treating scattering 
and bound-state problems. In the (renormalized) Numerov method, the bound 
states in a certain potential are calculated one at a time. To this end, we start 
with a trial energy and for this energy we try to construct a smooth wavefunc- 
tion that matches the boundary condition. If we fail, we need to have a criterion 
to improve our guess for the energy, so that a few times later we will succeed. 
We also need to be able to test if we have calculated all possible states, which is 
usually done by counting for each state the number of times that the wavefunc- 
tion goes through zero. However, it turns out that for large values of the order 
parameter, the functions u n i and v n e can have very different shapes from the 
noninteracting solutions, so that counting zeroes leads to problems and states 
can be missed. 

This problem is absent if we use a different approach. Namely, if we eval- 



uate the matrix H n £ of Eq. (1171 in a certain basis, then we can calculate 
all eigenstates and eigenenergies at once with matrix diagonalization. In Ref. 
|142j the harmonic oscillator basis was used, which analytically solves the di- 



agonal part of Eq. (115), while for the off-diagonal part integrals of the form 
(4>nem\A(r)\<f> n e m ) have to be calculated. We use another method that is based 
on the discrete variable representation (DVR). Using a DVR based on sine- 
functions (sinc-DVR) [234] , we have that the matrix H„£ is expressed in a basis 
that results in simple analytic expressions both for the kinetic energy operator 
and for position-dependent functions. 

We start with discussing the sinc-DVR method for the one-dimensional case, 
after which we generalize the method to the three-dimensional spherically sym- 
metric case. More details can be found in Ref. |234j . The basis set that is used 
is given by the following set of orthonormal functions 

Xj{x) = 7k sinc Ki~ J )]' (123) 

where sinc(a;) = sin(x)/x and j = —M, .., M with 2M + 1 the number of basis 
functions. These basis functions are to be evaluated on an equidistant grid 
Xi = iAx with i = —M, ..,M and grid spacing Ax. Note that Xj( x j) = 1/VAx, 
while for all other grid points Xj( x i) — 0- For the sinc-DVR basis, we have that 

Vtj = {Xi\V\Xi)tii, (124) 

for an arbitrary potential V(x), where the implied integral is to be evaluated 
using a quadrature rule on the discrete grid with constant quadrature weight 
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Wi = Ax. For the second derivative, we then have 



UX { Ax 2 (i-j) 2 ' 11 1 T J- 

As a result, we find that position-dependent potentials are diagonal in the sinc- 
DVR basis on the specified grid, while the kinetic energy operator has a simple 
analytic form. To illustrate the performance of the sinc-DVR method we may 
apply it to the one-dimensional harmonic oscillator. Using a step size of half 
a harmonic trap length, Ax = 0.5/, and using M — 10, we have to diagonalize 
a 21 x 21 hamiltonian matrix, for which standard numerical diagonalization 
routines can be used. Doing this, we get for the lowest twelve eigenvalues an 
error of less than 1 % compared to the analytic result E n = (n + 1/2)HQ. 

The three-dimensional harmonic oscillator in radial coordinates reduces for 
the spherically symmetric case to the corresponding one-dimensional problem 
after considering the functions 4>(r) = 4>{r)/r. For the range r € [0, oo], we can 
construct so-called wrapped sine- functions |234j . given by 

xf(r)=X j (r)±Xj(-r). (126) 

Here, we have j = 1,..,M and (ri,Wi) = (iAx,Ax) with i = 1,..,M for the 
radial grid and the quadrature weight. Note that the antisymmetric basis func- 
tions xj( r ) by construction fulfill the boundary condition xj(0) — 0, whereas 
the symmetric basis functions Xj{ r ) have a vanishing first derivative at r = 0. 
For the wrapped sinc-function basis sets, an r-dependent potential function has 
again only diagonal matrix elements, while the second derivative gives rise to 

r] 2 { 1 " 2 T 2 ( ~ 1)i+J ' if i-i 

/ Y ±|^_| Y ±\ = J 3 Ax' + Ax* (i+j)'> % ~3, (127] 

I Ax 2 (i-i) 2 + Ax 2 (i+j) 2 ' 11 / J' 

Applying the sinc-DVR method to the three-dimensional harmonic oscillator, 
we find that for I — the antisymmetric basis set x~ performs best, while for 
1=1 the symmetric basis set x + performs best. For higher angular momenta 
the difference between the two sets is very small. In the following, we choose 
the x~ basis for I even, and the x + basis for I odd. 

As a result, we are now in the position to evaluate the complete position- 



dependent 2x2 matrix H n < of Eq. (117) in the sinc-DVR basis. Using M grid 



points, we have for each angular momentum I a 2M x 2M matrix to diagonalize. 
Consider first the balanced case with h = 0. In the presence of the Cooper-pair 
condensate A(r), we obtain positive eigenenergies E n i for the quasiparticles 
with spin +, and the same negative energies — E n i for the quasiparticles with 
spin — . The opposite sign is because we reversed the order of the fermionic fields 



in Eq. (123 1. However, after the diagonalization we wish to bring the fermionic 
quasiparticles back in their normal order, so that then all quasiparticles are seen 
to have positive energies E n £. For the imbalanced case, we get after the normal 
ordering of the quasiparticlc fields that the dispersions are given by — ah + E n i 
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with a = ± the spin of the quasiparticle and E n i > 0. In the following, all 
summations are over positive energies E n g. 

For the gap profile A(r), we start with the profile given by the local-density 
approximation. To find a fully self-consistent solution to the BdG equations, 
the order parameter should satisfy 



A(r) 



4> + (t, r)0_(r,r)) 

E ^\ 1 u n i(r)v nl {r)[l - f(E ne + h) - f{E nl - h)] 

n,i 



(128) 



where in the second step we used Eq. (114 1 to write the expectation value 
in terms of the i/v fields. Since the action is diagonal in these quasiparticle 
fields, their correlation functions are readily calculated, resulting i n the Fermi 
occupation numbers /(e) = l/(e l3e + 1) in the second line of Eq. (128 1. The 
primed summation in Eq. ( |128[ ) indicates that we sum over positive energies, 
while The factor 21+1 arises from the summation over mi. Eq. (128) contains 



a divergence due to the use of the contact potential, as was discussed in Section 



2.5 In order to obtain a divergence-free expression |235j . we combine Eqs. (23 1 
and (fl28j) to 



(129) 



In the unitarity limit, we have that the left-hand side of Eq. (1291 is zero. 
The right-hand side of Eq. ( 129 ) is finite and we can find an analytic expression 



for the high-energy tail, which can be used to improve the convergence of the nu- 
merics. The analytic expression is proportional to A(r) and the proportionality 
factor T^" 1 (r) is calculated from 



A(r) = 
T A (r) 



-E«n(rR(r)0(£ n -£ A ) + l]r 



A(r) 



k 



l 

v 



E 

|k| = fc A (r) 

mA(r 



A(r) 



2(e k -M + ^ cx (r)) V 



A(r) 



2e k 



4ir 2 h 2 



2fc A (r)-fc M (r) log 



k A (r) + k^r) 



k A (r) ~ ^(r) 



(130) 



where we considered only high-energy states with energies above E Al since 9{x) 
is the Heaviside stepfunction, and where we introduced the wavevectors k^(r) = 
y / 2m((i - V cx (r))/h 2 and fc A (r) = y / 2m(E A - V cx (r))/h 2 . In the second step 



of Eq. (1301, we have used that for high-energy states, the gap profile A(r) 
and the potential ^ cxt (r) are slowly varying functions, so that we may use 
the semi-classical (i.e., WKB or local-density) approximation for the coherence 
factors M n (r) ~ Uk(r)/v / V and v n (r) Vk(r)/-\/V, where Uk(r) and Vk(r) are 
the analytically known homogeneous coherence factors evaluated with the local 
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Figure 29: The order-parameter profile for a spin-polarized superfluid in the unitarity limit 
calculated with the Bogoliubov-de Gennes method at zero temperature. The figure was pro- 
duced using fi = 32hH! and h = 2Ahui, which results in a total number of particles given 
by TV = 3 X 10 4 and a total polarization of P = 0.95. The horizontal axis was scaled by 
R+ = \/2/i+/ma} 2 . The interface between the normal and superfluid phase is seen to have a 
nonzero width, while near the interface oscillations in the order parameter are observed. 



chemical potential fi — V cx (r), the local order parameter A(r), and the local 
wave vector k = ^/2m(E n ~ V cx (r))/h 2 . These homo geneous coherence factors 
were given below Eq. (|28| and for high energies they reduce to 



2 Uk (r> k (r) = AW » AW (131) 



which was used in Eq. (130). By combining Eqs. (1291 and (130) we finally 



obtain the appropriate gap equation in the unitarity limit 

W«n(r)<(r) 

1 A ( T ) 



(132) 



which, together with Eq. (115), gives the mean-field BdG approach. The self- 



consistent solution for the imbalanced Fermi gas is presented in the next section. 



6.2.3. Results 
In Fig 



29 



we show the results for the order parameter profile by applying 
the Bogoliubov-de Gennes equations to the imbalanced Fermi gas. We used 
(j, = 32/kI>, h = 24to, E A = 50/i, Ax = 0.05/, and M = 250. The local-density 
approximation predicts for these chemcial potentials a sharp interface at a radius 
of about r = 21. Around this radius, the BdG approach is indeed seen to give a 
smooth interface, where the most remarkable feature is the oscillatory behavior 
of the order-parameter close to the interface. This oscillatory behavior of the 
gap has been studied for example near superconducting-ferromagnetic interfaces 
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and is known in this case as the proximity effect. It follows from a selfconsistent 
study of Andreev reflection from the interface [2361 12311 1152] , where reflected 
holes have a different wavevector than incoming particles due to a difference 
in chemical potentials. Although the proximity effect shares similarities to LO 
type pairing, since both involve pairing correlations at nonzero center-of-mass 
wavevectors, we recall that the LO phase requires the spontaneous formation of 
a oscillatory superfluid, as happens below a Lifshitz point. This is in contrast 
to the present case, where the oscillations occur due to the presence of a sharp 
interface, as happens below a tricritical point. 

In order to study the resulting density profiles, we use the following equations 
for the densities 



n±(r) = (01(T,r)0±(T,r)) 

'i£(r)[l - f(E n ±h]+ u 2 n (r)f(E n T h), 



(133) 



which are similar to Eq. (28) for the homogeneous case. The high-energy tail of 



Eq. (1331 can be simplified to 



Y^vl{v)9{E n - E A ) 



£ 

|k|=fc A (r 



A(lf 



4(e k -^ + F oxt (r)) 2 ^ 



(134) 



with the same cutoff A as for the gap equation. The resulting integral on the 
right-hand side can be performed analytically. The chemical potentials lead to 
a total number of particles of N = 3 x 10 4 and a total polarization in the trap of 
P 



0.95. The resulting density profiles are shown in Fig. 30 



We see that the 

densities also show the oscillatory behavior close to the interface. So far, we have 
ignored the inclusion of self-energy effects, so that we do not have the correct 
equations of state for the normal state and superfluid state. As a result, the 
interface is not in the correct position compared to experiments. In principle, 
we could include the information about the interaction effects from Monte-Carlo 
calculations in a similar way as for the Landau- Ginzburg approach, but then 
the densities and the gap would have to be determined by the appropriate 
derivatives to the thermodynamic potential, so that Eqs. ( 132 ) and ( 133 1 would 



not be valid anymore. Such a generalized BdG approach was put forward in Ref. 
|237j . where also a DVR-based method was used to solve the BdG equations. 
Since the prediction for the proximity effect is so far only based on the mean- 
field BdG equations, the existence of the effect in the unitarity limit remains 
uncertain until confirmed or invalidated by experiments. 



7. Discussion and outlook 

In this review, we have treated the imbalanced atomic Fermi gas with strong 
attractive interactions, for which pairing plays an important role. We have fo- 
cussed mainly on the population-imbalanced case, for which in recent years 
many exciting experiments and a large amount of theoretical calculations have 
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r/R+ 

Figure 30: The majority and minority density profiles of the polarized superfluid in the 
unitarity limit at zero temperature calculated using the Bogoliubov-de Gennes method. The 
same parameters were used as for Fig. |29| The vertical axis is scaled using no = n+(0). 
Expecially in the majority density profile (the upper curve) clear oscillations are seen near 
the interface. 



been preformed. The large amount of interest in the topic started in the be- 
ginning of 2006, when two experimental groups obtained a full control over the 
spin imbalance or polarization P in an ultracold atomic Fermi gas [241 125) . This 
made it possible to study experimentally in detail the fundamental question: 
what happens to a superfluid Fermi mixture upon frustrating the pairing by 
increasing the spin polarization? We have seen that in the atomic Fermi gas 
there are only two-body attractions between particles of opposite spin, so that 
by increasing the polarization, also the number of particles that cannot pair 
increases. As a result, we argued that there must be a transition from the 
superfluid to the normal state somewhere between zero and full spin polariza- 
tion. But when? And how? One of the experimental and theoretical findings 
was that these questions depend on the temperature of the Fermi gas. Namely, 
when the temperature is above the tricritical point, the transition is smooth 
(or second-order) and occurs at lower polarizations, while below the tricritical 
point, the transition is sharp (or first-order) and occurs at higher polarizations. 
We have seen that this behavior can be qualitatively explained by mean-field 
theory, where we mapped out the phase diagram both for the homogeneous and 
the trapped case. However, on a quantitative level the mean-field critical polar- 
ization, or Chandrasekhar-Clogston limit, turned out to be very different from 
the experimental results. 

This problem was solved by Monte-Carlo studies of the strongly-interacting 
normal state and the superfluid state at unitarity, which resulted also in very 
good quantitative agreement with experiments [150 . Although in this review we 
have not reviewed Monte-Carlo techniques, we have discussed several other ways 
to go beyond mean-field theory, such as including Gaussian order-parameter 
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fluctuations (the random-phase approximation) , diagrammatic many-body scat- 
tering calculations, and renormalization group calculations. Although these 
methods are in principle less accurate than Monte-Carlo simulations, they have 
also some advantages. First, they require much less numerical effort. Second, 
they show more clearly the relevant physics involved. And third, they can be 
used at nonzero temperatures and polarizations, while Monte-Carlo simulations 
have so far been restricted to either zero temperature or to the balanced case. 
Together, the theoretical approaches allow for a rather complete present under- 
standing of the equilibrium properties of strongly-interacting imbalanced Fermi 
gases. This holds both for the homogeneous case, and for the trapped case with 
the use of the local-density approximation. 

In this review, we also calculated the mean-field phase diagram for the mass- 
imbalanced case of the experimentally relevant 6 Li- 40 K mixture, which not only 
contains a tricritical point in the unitarity limit, but also a Lifshitz point, so 
that supersolidity is expected to occur. Key open questions in this respect 
are the precise phase structure below the Lifshitz point and a determination 
of the crystalline structure of the various supersolid phases that are realized. 
Moreover, the Sarma regime occupies a large part of the phase diagram for 
the 6 Li- 40 K mixture, so that experimental observation of a gapless superfluid 
with (smoothened) Fermi surfaces might be within reach. Finally, we treated 
the trapped inhomogeneous Fermi gas beyond the local-density approximation 
in order to study in more detail the supcrfluid-normal interface that occurs 
below the tricritical point. With a Landau-Ginzburg approach we calculated the 
surface tension of the interface, while with the Bogoliubov-de Gennes approach 
also order-parameter oscillations near the interface were observed due to the 
proximity effect. However, in order to compare these theoretical results for the 
interface in detail with experiments more accurate experimental data for the 
interface would be required. 

Although significant steps have been made in understanding imbalanced 
Fermi gases, there is still room for more research. From the theoretical side, 
further improvements could be made by having also Monte-Carlo data for the 
imbalanced Fermi gas at nonzero temperatures, and in particular for the tri- 
critical point. Another theoretical improvement of interest would be a more 
detailed account of screening effects of the interaction by particle-hole excita- 
tions, since screening is usually either neglected in analytical approaches, or 
treated in simple approximations. The latter typically ignore the precise mo- 
mentum and frequency dependence of the particle-hole excitations. From the 
experimental side, a more careful study of the superfluid-normal interface is 
desirable. When the local-density approximation applies, this would then allow 
for a more careful experimental mapping of the homogeneous phase diagram, 
which up to now contains very few points with large error bars p. 48] . Seen the 
experimental accuracy that is currently being reached for balanced Fermi gases 
[711 1173] , more precise data for the phase diagram of the imbalanced Fermi gas 
should be within reach. This would then also result in a more precise determi- 
nation of the tricritical point. If the interface is studied more accurately, then 
it is also natural to check if the proximity effect can be observed in a cold-atom 
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system. 

Furthermore, an experimental determination of the surface tension would be 
a test for our knowledge of the superfluid state beyond equilibrium, since the 
surface tension probes the unstable barrier of the superfluid grand potential as 
a function of the order parameter. It would also be interesting to study vari- 
ous other non-equilibrium aspects for imbalanced Fermi gases. These systems 
contain a first-order phase transition, so that they are ideally suited to study 
nucleation in a quantum system. Novel topics of interest for imbalanced Fermi 
systems include transport properties, artificial gauge fields, spin-orbit coupling, 
optical lattices, p-wave interactions, reduced dimensions, and so on. These top- 
ics are all expected to become experimentally available for imbalanced Fermi 
gases, or have already recently been realized, so that imbalanced Fermi gases 
are expected to keep us interested for many years to come. 
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